*1}v*K_vu~"=d?2+1u˺zN?KX4Ur(|v'gYt:EIX__yˉ*=ϳ\b9 mޟw9 &
O >
C"e@* VӐ2 * E9kҨQ$QGD/2ICfC aпIS7܊tjV#l c 4Ng#HzmM*ꌸG$U!T@\]6h5Ԑ81Hh=.h44*b)5sǗrP%Coج6"z[Dyp+zs;!Q[ (?JiX!mV#:H
<hMdT/sոMNMž/h$xXڈSg_RfU\22Д.'\p|;`X_Po҆M:w˄RM
jr.p ϽH.5A5`'L=(8&q02C, V9wX5܀1-ldRXM@v *t1 < vJG% q/0lOAZIeo1 W*L M>d<K 6Z MoS-wM( ۹P mݿ,
N^J)f:*4]p^T720Q%L% yi+h`,+zUtzR|8\z |˘>( ീ~S0dqȘ/TfNjMAn 2f6Yiϯds2f<_Y_
pCKg;qα>1P}t"tBCi5#Ih\o6t~Kd?kr};Čf4U.??[j}Jj*ւ!,JG飤o VfY֑++
aoζLKGkvDy\
e&U4;[UOy 77Yr z1ZE
Ȉnkn:7 eJ E}Hszc N$F}Hz:^Qݴ,F4?ݞ1z mmFVuIe!ݦ
5Q4v}ZfLh%ytJ:5_a b# {1ۆVpD*H$bZ[fh6c9kuVժ\>/yB;hSP4tRА%@ C4>"S-B
k>L$TJsoúθS\JՕ
wqq.Sʺ;V׆9oJ0 sNa֊+i)kDk~^Q5 T^[voDWPZK}!r b 37K34lLǗï)CuҼhjz#co^< %kez
5 DBKBEV2DR[y-2{jނfXn ef|
ŋ4|X7_/h!` 1D =7:čƌ=M]YGA;34VLk5.ā:\tH_(OQ #6Gpz?$5T8^K]zS9Uaי/7<8 7<Şa:< ;EOV%'WQ} 1wr6Yb29M; ֨9K`*sDk5/CaRT-x&b,VGI6'`m|N F)aH-w!) Ml̶|<&|5uئ+KnuM
$!L,l;TzUl9 {7]UUհsä $U-PÀ߄u@sO x!Fɠ˚Sq(ֻt j5uIX1nOYŠ|;qZZ[\
x6|f Or
ϵt2hx[>ۚS`a9LiH s8ڤ^ k?m "z
r5m$AKBipG)"w7isTa:Âz+kM(b!Z]i5FicԥTVj&)OMJai(äaϨ3E| Hѹ 'LRf~<Xn%iA? J8ݳuW[K /-A
i3q& f7\{rqZl(!ur>5,Ɔ@ jGE-[1/@۱yH+%6ln.ҥT$ %cN-CZ 2/'L,7N!K,v/Fg 9 8]1-n ҬFWhѾWA`f.52q* \ x\D'CliٵMtmFMx:J㨣\.nܙ_KGO(w4r c8nĨ;׃3b5^ʿ>a !Wn pG4_ v(wPc>&3֠1A(Oi)9'' y&7H I[y12v"e̤IpXG,83RE HJPљ 8c挲D5` "- idԉ?jHS& &8H "pz[MЩxEM,-iX%c
@mmΚ +ۤt !U]uҭŁMU7gT#qE t{U3\U@ҵ)g@˘rJQ"ޕ~C(: z&py uy 7 6;*rSw D6A3tPJE<:dߦ7@|uji$EiXBB,^IXO%>tdY8i!n .&,~ ۜ tt>L#1mi
U fQ)},UwF$ZjL F%#tj䐴OKU}W-qM]'Txa+Ъ>Wby:$ԜǓh[ I8 w'`mVr<ݿ+gWM2~"xA @2l4A^"Fz3ٟߚ &//*SyT ^mqjt3[qJꖇmSΒ =7L7̜PzuaG?GXULQIE0`<=V.zyF&@ZxVrxQX7\4jFKLۆ% ONKr̫C:STrl$s<[u}y"[({6Yy 7 :
M5u X'9JiӺ)rNCgzXCYgy/9di=JYr '?Ku.m|nϒ6moQ
9|H 9d'VE#<_$ErxY$K#3_mBa;wb83*mz7)D.v{"Paۂ3ڍTH%Mnv{mM]SJ]즻Y|%;}n[ق:ߖ^#ݕnxtz=mts@:ގ!h`Wkh7 PK 8BF
FuncDesigner/sle.pyWKo6W v l((`Sb-mbiR%)ޢ3eMtp$frXoɕF{˱%3CĆϘepx{cVHkN]m}é@ksPқ^^\J8NO{,-AIS 3^
%`Sa mr3jh#QyPaG0K2'D][Ccv.۽%ғ$W~,YiM44`K Ak-uZ*j"b R+;U.\ zӫ>Cޤhw+ױ$̜9sw3lwe}>M?eMTQOBߡġAM5h.\J@izmJ;FZ' E\U,TY[0㳾 8Ϻ[KRl?ʳan:0? 3Vyi䌲IXt
[MvpE_Z=j&C"7aeR]U >M PI[dCiA T<5YoDK
X N )-u~ϫ:0bҠ ol#< 4qΑ+&seYTΕzt hU>Zoahj_H6n,Ft$l(إ b`gKlŰF{aH (U3YT}^yWI h gLBT!)1(+ } GdyejSi zyU |Hg#lW$:5_g6ř7L22E:k8I`t_sDG'A_M ]*p;1۟٤82+>?ʑM<"ː t% rsPIu6Ӱgwｿ]^-= nz]9K:>H
?^̧<3 i Ii ^m r z`#wLt? XKgO9n/Af =bȵ IIbAyF&~t!W\q7wbhF ذߊe @:8oU4Ԇj(W9kj]7@5_=+5WojC /!)04TwH;`f/>o&M,QD'cbPo00
=8)h C"%e3g1`]x1*=tF$̥}hw Yv@^_)3dh QuG>RN@DT.QYw9 UX{ 0:-}ێ|9S֗|`d&aRi {GoXU6 :+xYeڻ!@]\8TlB
:Y2 kPPZ mݏ$.ݙ:ؠbBCjpt YqA%4i`]ILb>ܢO0C:?V& :
2.';BKs͖`3= ?}bt6y[fC-Ti@: Lx#qG b:xG~g ߀,z Ā h C敩,xtP $)%dfh]S#3/"_&tW>5T mf. hR E̱~iD.^i/ _]dX I=AD:Vp rFcE`Xs0͝\Ug"z^Zg ˞Mw=&?%~~L#Bgǫoߑx D'o0ƇCWBNX{H!'3sLudT#LL'xl db >p\oa {n6APK 8B
X FuncDesigner/stencils.pycUKo#E q.0/HXBr Eb'f4dW߁\9# ؙUJ;ߍ~ ~}B(Q4P84pHAP\ĥk%=J4(٥Ҟltݭp tqG^b1cd]eY> M̹" =h^Y Mb47;cMNT(zhY¡ 37-]*\MTh)1LW=RܵwTo GWiBp%Pm#ZL SA+/)%Ci=j-+/WQYm NU=*8F*ݖHאA-Tzh=GihG֩48}d- Y_XԤ!hs"UǒؐݫNYNNDoQQ[%'juџW9:;t9Bsf,ns__4TO醥4=g-oĞwa ͎cgџP=Fc*#bR@ |a2۰< b py cQv qm %!Yfybq7>cqbݻQ_$4{}~'W#qosJlٝ`H5pj|ԡsw*_]S|;U/X.22 E9(adfdA2Axv`Pqe+K@Fi3xvjgX4\Al@ 3 3.LOZZzYFE]]QEu]|ݺUw+h@vu"؇88\&}ڦ9iKiQ~V=9HnKW<ufH^XjUObYt
f
FHtBD(}%%qv彗qxȞy2dņy|#* JU9wtz.En_u#~ΦMk!zXڈG8|qX--O+1)*ϛ;*н3<%X7TnHDҍP"6^RH73~Q!{B%;ZD>Nmsк!+R3b@KE\4
Cx&4{8civHjx|dވ͍3*tL-Lڕh,DFM'Yl
DAOy@ _}:-~ -"U1۠~_o ;XYs}j*Ljt0&-F-
1.7
v +k h4
6ԒZOegFV,9)E=K ?|?:ei/)@ɲwO4lV]w}{-: 3|SaDM|w)K7Z0G3teVZ+V{@0K &kdqaS"4 x=| d FM vW֯ yEL]#oF mEak]){j$y`'R;B&ĭ'QAI~ 8r#q^ũ!'H\8
VYE
*Vi lY$FlLjz{
DL#F( j I wΠ DAYԲPH@;kV8 $7ztbunonb^c88@֦yp ;}0W)d.&ޣ'6 㱘c ^n< 0w: J709ž?\W7
66p ^RwCoH8t);}:+1y>Ŧ/ 09\N륧_0O&cw^sCu=ԏbl,;:/Df5$A.cQM(?637[~i &'!!x{KǡXsU*M @S X87Y{00> )bIA
0W{wEhO1UD犥;n4A~Usev^6_/kt|:2eN67oXĖnX =4aj{ WNj*Q^ b6žˈ}ثPK 8BurD
: FuncDesigner/translator.pycXs G=Jϕ,m8k;-[@cNJYW[eӒFڙYʒkb'8 J_ŕ*rKs@q 띕Pp~u956
\3 S!5E5 \jYhcYDE L2EԜhUYhSMSsjP4KY9G
STf>Yh \E
%
Ʃ2ڇ⟕91 t5ȳ=8OlOS 1ᛎ3ӺyeȔR#&Ɓ]-&LXg""̲S# h=#kƊA9ʜz#Qd@q58! ?t ?>v^cE =s]EO JKyPֶl{c #v ~I[-~ЃUS>T(?TŘWL;]+R
=VHOa =!nfzb ce(_P />ѕ
rS#$&3df1e9 8$Dž^X T6_HK5T,%AiԧL[,նFm"
9ufrs^-OCιS@r/@^ "ǤdaU"nk+;3w:'dCA7iONo{̛'ASJ܋t-L[~;LMV`w̴(u^O6H.{e?q_tutgK$̹Ɓ> ez\j_i" LeUTꔚ,Δ*al_La>rL a V!XV& _ 7=FJF54MRq .rj{_ٗb#2
{q aŻ2=` R?H?-* /7 xqxrd/ y_ \wD̈ikby{3lzpgղέ⬚K I^b3i6!!p 4r2a=iu#HƂI7oQ@:Jvɾt{/
>pe{eJٽk 0_m#X0\adJß!yYQ[vKu,)"薈h ~/ `cêkÂ8\GN\循x;CZU O
/5't<
1Ja~YJ
>Xu&r)R|W.6cHO`dL
L8MA;Y*џЀ5LўK_Uvrd}s,3
,Si15+d^.2OvOq "P25xzk5/Sjo ^H;UkP;ǜs^ yüoOvx~C
9 yjH0 ~U/^pt`K~v~C#;: `YI PVpsUg3w/A7y`x>V -.|qG^>o)?ipAJk ¸)]=U~Enpn$ w8z墣oҰaҫg;I0>o$Ad{Xс49 mZר1Zv'RS_
BW$a=7(|䣆GtK -4,7 @iΣڬ %6!6ϟ߲o뒛J5o!zx6Bac{o% շMX m?ly_˫EhtyCO[n9qu(A[\m@rLf]1ށwJy^!xxՔӷ+
)q{4$7|4 V,jL^6uXZ:^-0"u !Z.}Bbmj9/3baɂu
Ր;"-iqL&p`Eۅk}[5:ޫPú<:Y[nUU^ufWhf3q.ӔEgoMUzY >Q:gPEN8jr jb?ҭmx[6[-ruxญa&: D|m,eVU/q/L5ۏI{*]Xa9m&q%6कFˋ}±[Qm{ 7__YЯ;캳j PK 8B2 EGG-INFO/dependency_links.txtPK \8B2 > EGG-INFO/not-zip-safePK 8BVg o t EGG-INFO/PKG-INFOPK 8BSv
EGG-INFO/requires.txtPK 8BNƇ J EGG-INFO/SOURCES.txtPK 8BaD<
P EGG-INFO/top_level.txtPK 8BWP+ 0
FuncDesigner/__init__.pyPK 8Bb_;
FuncDesigner/__init__.pycPK 8Bu / f FuncDesigner/baseClasses.pyPK 8Bb FuncDesigner/baseClasses.pycPK 8Bzے2 S FuncDesigner/dae.pyPK 8Ba1 FuncDesigner/dae.pycPK 8Bh 9 FuncDesigner/FDmisc.pyPK 8B E ! % FuncDesigner/FDmisc.pycPK 8Bl X
2 FuncDesigner/integrate.pyPK 8BGB
?7 FuncDesigner/integrate.pycPK 8B 8 S< FuncDesigner/interpolate.pyPK 8Bo
H i SB FuncDesigner/interpolate.pycPK 8B0 1:
return ooarray([IMPLICATION(condition, elem) for elem in args]) if condition is not True else args
return NOT(condition & NOT(args[0])) if condition is not True else args[0]
ifThen = IMPLICATION
from sle import sle
from ode import ode
from dae import dae
from overloads import *
from stencils import d, d2
#from overloads import _sum as sum
from interpolate import scipy_InterpolatedUnivariateSpline as interpolator
from integrate import integrator
isE = False
try:
import enthought
isE = True
except ImportError:
pass
try:
import envisage
import mayavi
isE = True
except ImportError:
pass
try:
import xy
isE = False
except ImportError:
pass
if isE:
s = """
Seems like you are using OpenOpt from
commercial Enthought Python Distribution;
consider using free GPL-licensed alternatives
PythonXY (http://www.pythonxy.com) or
Sage (http://sagemath.org) instead.
"""
print(s)
del(isE, curr_dir, os, sys)
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/baseClasses.py 000644 000766 000024 00000001457 12100431107 026216 0 ustar 00rfadams staff 000000 000000 from numpy import ndarray
class OOArray(ndarray):
pass
class MultiArray(ndarray):
pass
def distrib_err_fcn(*args, **kw):
from FDmisc import FuncDesignerException
raise FuncDesignerException('''
direct operations (like +, -, *, /, ** etc) on stochastic distributions are forbidden,
you should declare FuncDesigner variables, define function(s) on them
and then get new distribution via evaluating the obtained oofun(s) on a data point
''')
stochasticDistribution = 'stochastic distribution'
class Stochastic:
#__array_priority__ = 100
def __init__(self):
self._str = stochasticDistribution
__repr__ = lambda self: self._str
__add__ = __mul__ = __pow__ = __rpow__ = __rmul__ = __radd__ = __neg__ = __pos__ = distrib_err_fcn
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/dae.py 000644 000766 000024 00000004437 12100431107 024520 0 ustar 00rfadams staff 000000 000000 from ooSystem import ooSystem as oosystem
from FDmisc import FuncDesignerException
from numpy import ndarray
class dae:
def __init__(self, equations, time, startPoint = None):#, **kw):
self.equations = equations
self.startPoint = startPoint
self.time = time
s = 'for DAE time must be dict of len 1 or array, '
if type(time) == dict:
if len(time) != 1:
raise FuncDesignerException(s + 'got dict of len ' + str(len(time)))
self.timeInterval = asarray(next(iter(time.values())))
self.N = self.timeInterval.size
else:
if type(time) not in (list, tuple, ndarray):
raise FuncDesignerException(s + 'got type %s insead ' + str(type(time)))
self.N = len(time)
self.timeInterval = time
if self.N < 2:
raise FuncDesignerException('lenght of time must be at least 2')
# !!!!!!!!!!
# TODO: freeVars, fixedVars
# !!!!!!!!!!
# for fn in ['freeVars', 'fixedVars']:
# if fn in kw
def solve(self, solver = None, **kw):
S = oosystem()
S &= self.equations
if self.startPoint is not None:
startPoint = self.startPoint
else:
Dep = set()
Dep.update(*[eq._getDep() for eq in self.equations])
if None in Dep:
Dep.remove(None)
startPoint = {}
for v in Dep:
if 'size' in v.__dict__:
startPoint[v] = v.size
else:
startPoint[v] = [0]*self.N
#startPoint = dict([(v, [0]*v.size) for v in Dep])
kw2 = kw.copy()
if solver is not None:
kw2['solver'] = solver
r = S.solve(startPoint, **kw2)#, iprint=1)
r.plot = self.plot
self.r = r
return r
def plot(self, v, grid='on'):
try:
from pylab import plot, grid as Grid, show, legend
except ImportError:
raise FuncDesignerException('to plot DAE results you should have matplotlib installed')
f, = plot(self.timeInterval, self.r(v))
legend([f], [v.name])
Grid(grid)
show()
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/FDmisc.py 000644 000766 000024 00000017471 12100431107 025136 0 ustar 00rfadams staff 000000 000000 from numpy import asscalar, isscalar, asfarray, ndarray, prod
import numpy as np
from baseClasses import MultiArray
scipyInstalled = True
try:
import scipy
import scipy.sparse as SP
except:
scipyInstalled = False
from baseClasses import Stochastic
class FuncDesignerException(BaseException):
def __init__(self, msg):
self.msg = msg
def __str__(self):
return self.msg
def checkSizes(a, b):
if a.size != 1 and b.size != 1 and a.size != b.size:
raise FuncDesignerException('operation of oovar/oofun ' + a.name + \
' and object with inappropriate size:' + str(a.size) + ' vs ' + b.size)
scipyAbsentMsg = 'Probably scipy installation could speed up running the code involved'
pwSet = set()
def pWarn(msg):
if msg in pwSet: return
pwSet.add(msg)
print('FuncDesigner warning: ' + msg)
class diagonal:
isOnes = False
__array_priority__ = 150000# set it greater than 1 to prevent invoking numpy array __mul__ etc
def __init__(self, arr, scalarMultiplier=1.0, size=0):
#assert arr is None or arr.ndim <= 1
self.diag = arr.copy() if arr is not None else None # may be None, then n has to be provided
self.scalarMultiplier = scalarMultiplier if isscalar(scalarMultiplier) \
else asscalar(scalarMultiplier) if type(scalarMultiplier) == ndarray\
else scalarMultiplier[0, 0] if scipyInstalled and SP.isspmatrix(scalarMultiplier)\
else raise_except()
self.size = arr.size if size == 0 else size
if arr is None:
self.isOnes = True
copy = lambda self: diagonal(self.diag, scalarMultiplier = self.scalarMultiplier, size = self.size)
def toarray(self):
if self.isOnes:
tmp = np.empty(self.size)
# for PyPy compatibility
scalarMultiplier = asscalar(self.scalarMultiplier) if type(self.scalarMultiplier) == ndarray else self.scalarMultiplier
tmp.fill(scalarMultiplier)
return np.diag(tmp)
else:
return np.diag(self.diag * self.scalarMultiplier)
def resolve(self, useSparse):
if useSparse in (True, 'auto') and scipyInstalled and self.size > 50:
if self.isOnes:
tmp = np.empty(self.size)
tmp.fill(self.scalarMultiplier)
else:
tmp = self.diag*self.scalarMultiplier
return SP.dia_matrix((tmp,0), shape=(self.size,self.size))
else:
return self.toarray()
def __add__(self, item):
if type(item) == DiagonalType:
# TODO: mb use other.diag.copy(), self.diag.copy() for more safety, especially for parallel computations?
if self.isOnes and item.isOnes:
return diagonal(None, self.scalarMultiplier + item.scalarMultiplier, size=self.size)
else:
if self.isOnes:
d1 = np.empty(self.size)
d1.fill(self.scalarMultiplier )
else:
d1 = self.diag
if item.isOnes:
d2 = np.empty(item.size)
d2.fill(item.scalarMultiplier )
else:
d2 = item.diag
return diagonal(d1 * self.scalarMultiplier + d2 * item.scalarMultiplier)
elif np.isscalar(item) or type(item) == np.ndarray:
return self.resolve(False)+item
else: # sparse matrix
assert SP.isspmatrix(item)
return self.resolve(True)+item
def __radd__(self, item):
return self.__add__(item)
def __neg__(self):
return diagonal(self.diag, -self.scalarMultiplier, size=self.size)
def __mul__(self, item):
#!!! PERFORMS MATRIX MULTIPLICATION!!!
if np.isscalar(item):
return diagonal(self.diag, item*self.scalarMultiplier, size=self.size)
if type(item) == DiagonalType:#diagonal:
scalarMultiplier = item.scalarMultiplier * self.scalarMultiplier
if self.isOnes:
diag = item.diag
elif item.isOnes:
diag = self.diag
else:
diag = self.diag * item.diag
return diagonal(diag, scalarMultiplier, size=self.size)
elif isinstance(item, np.ndarray):
if item.size == 1:
return diagonal(self.diag, scalarMultiplier = np.asscalar(item)*self.scalarMultiplier, size=self.size)
elif min(item.shape) == 1:
#TODO: assert item.ndim <= 2
r = self.scalarMultiplier*item.flatten()
if self.diag is not None: r *= self.diag
return r.reshape(item.shape)
else:
# new; TODO: improve it
if self.isOnes:
D = np.empty(self.size)
D.fill(self.scalarMultiplier)
else:
D = self.scalarMultiplier * self.diag if self.scalarMultiplier != 1.0 else self.diag
return D.reshape(-1, 1) * item # ! different shapes !
# T = np.dot(self.resolve(False), item)
# from numpy import array_equal, all
# assert array_equal(T.shape, T2.shape) and all(T==T2)
# print '!'
#prev
# !!!!!!!!!! TODO: rework it!!!!!!!!!!!
# if self.size < 100 or not scipyInstalled:
# return np.dot(self.resolve(False), item)
# else:
# return self.resolve(True)._mul_sparse_matrix(item)
else:
#assert SP.isspmatrix(item)
if prod(item.shape) == 1:
return diagonal(self.diag, scalarMultiplier = self.scalarMultiplier*item[0, 0], size=self.size)
else:
tmp = self.resolve(True)
if not SP.isspmatrix(tmp): # currently lil_matrix and K^ works very slow on sparse matrices
tmp = SP.lil_matrix(tmp) # r.resolve(True) can yield dense ndarray
return tmp._mul_sparse_matrix(item)
#return r
def __getattr__(self, attr):
if attr == 'T': return self # TODO: mb using copy will be more safe
elif attr == 'shape': return self.size, self.size
elif attr == 'ndim': return 2
raise AttributeError('you are trying to obtain incorrect attribute "%s" for FuncDesigner diagonal' %attr)
def __rmul__(self, item):
return self.__mul__(item) if isscalar(item) else self.__mul__(item.T).T
def __div__(self, other):
#TODO: check it
if isinstance(other, np.ndarray) and other.size == 1: other = np.asscalar(other)
if np.isscalar(other) or prod(other.shape)==1:
return diagonal(self.diag, self.scalarMultiplier/other, size=self.size)
else:
# TODO: check it
return diagonal(self.diag/other if self.diag is not None else 1.0/other, self.scalarMultiplier, size=self.size)
DiagonalType = type(diagonal(np.array([0, 0])))
Eye = lambda n: 1.0 if n == 1 else diagonal(None, size=n)
def Diag(x, *args, **kw):
if isscalar(x) or (type(x)==ndarray and x.size == 1) or isinstance(x, (Stochastic, MultiArray)):
return x
else:
return diagonal(asfarray(x) if x is not None else x, *args, **kw)
class fixedVarsScheduleID:
fixedVarsScheduleID = 0
def _getDiffVarsID(*args):
fixedVarsScheduleID.fixedVarsScheduleID += 1
return fixedVarsScheduleID.fixedVarsScheduleID
DiffVarsID = fixedVarsScheduleID()
_getDiffVarsID = lambda *args: DiffVarsID._getDiffVarsID(*args)
try:
import numpypy
isPyPy = True
except ImportError:
isPyPy = False
def raise_except(*args, **kwargs):
raise FuncDesignerException('bug in FuncDesigner engine, inform developers')
class Extras:
pass
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/integrate.py 000644 000766 000024 00000006530 12100431107 025745 0 ustar 00rfadams staff 000000 000000 from ooFun import oofun, atleast_oofun
import numpy as np
from FDmisc import FuncDesignerException
from translator import FuncDesignerTranslator
from ooPoint import ooPoint as oopoint
from ooVar import oovar
try:
from scipy import integrate
scipyInstalled = True
except:
scipyInstalled = False
def integrator(func, domain, **kwargs):
if not scipyInstalled:
raise FuncDesignerException('to use scipy_integrate_quad you should have scipy installed, see scipy.org')
# if not isinstance(domain, dict) or len(domain) != 1:
# raise FuncDesignerException('currently integration domain must be of type Pythoon dict {integr_var:(val_from, val_to)} only')
#integration_var = domain.keys()[0]
#_from, _to = domain.values()[0]
integration_var, a, b = domain
if not isinstance(integration_var, oovar):
raise FuncDesignerException('integration variable must be FuncDesigner oovar')
a, b, func = atleast_oofun(a), atleast_oofun(b), atleast_oofun(func)
def f(point=None):
p2 = point.copy()
#if integration_var not in p2:
#p2[integration_var] = 0.0 # to ajust size, currently only R^1 is implemented
#T = FuncDesignerTranslator(p2)
def vect_func(x):
p2[integration_var] = x
tmp = func(p2)
if np.isscalar(tmp): return tmp
elif tmp.size == 1: return np.asscalar(tmp)
else: FuncDesignerException('incorrect data type, probably bug in uncDesigner kernel')
#vect_func = lambda x: func(T.vector2point(x))
# TODO: better handling of fixed variables
return integrate.quad(vect_func, a(point), b(point), **kwargs)[0]
# def aux_f(point):
# point = oopoint(_point)
# if a(point).size != b(point).size: raise FuncDesignerException('sizes of point-from and point-to must be equal')
# if a(point).size != 1 or b(point).size != 1: raise FuncDesignerException('currently integration is implemented for single variable only')
# dep = a._getDep() | b._getDep() | func._getDep()
# involved_vars = set(point.keys()).add(integration_var)
# if not dep.issubset(involved_vars): raise FuncDesignerException('user-provided point for integration has no information on some variables')
#
# p2 = point.copy()
# if integration_var not in p2:
# p2[integration_var] = 0.0 # to ajust size, currently only R^1 is implemented
# T = FuncDesignerTranslator(p2)
# vect_func = lambda x: func(T.vector2point(x))
# TODO: derivatives
# !!!!!!!!!!!!! TODO: derivatives should not be zeros! Fix it!
r = oofun(f, None)
r.fun = lambda *args: f(point=r._Point)
# TODO : use decorators here
tmp_f = r._getFunc
tmp_D = r._D
def aux_f(*args, **kwargs):
if isinstance(args[0], dict):
r._Point = args[0]
return tmp_f(*args, **kwargs)
def aux_D(*args, **kwargs):
raise FuncDesignerException('derivatives from scipy_quad are not implemented yet')
if isinstance(args[0], dict):
r._Point = args[0]
return tmp_D(*args, **kwargs)
r._getFunc = aux_f
r._D = aux_D
return r
# r = oofun(lambda x, y: integrate.quad(func, x, y, **kwargs)[0], input = [a, b])
# r.d = (lambda x, y: -func(x), lambda x, y: func(y))
# return r
#def
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/interpolate.py 000644 000766 000024 00000011070 12100431107 026304 0 ustar 00rfadams staff 000000 000000 from ooFun import oofun
import numpy as np
from FDmisc import FuncDesignerException
try:
from scipy import interpolate
scipyInstalled = True
except:
scipyInstalled = False
def scipy_InterpolatedUnivariateSpline(*args, **kwargs):
if not scipyInstalled:
raise FuncDesignerException('to use scipy_InterpolatedUnivariateSpline you should have scipy installed, see scipy.org')
assert len(args)>1
assert not isinstance(args[0], oofun) and not isinstance(args[1], oofun), \
'init scipy splines from oovar/oofun content is not implemented yet'
S = interpolate.InterpolatedUnivariateSpline(*args, **kwargs)
return SplineGenerator(S, *args, **kwargs)
# TODO: check does isCostly = True better than False for small-scale, medium-scale, large-scale
# return SplineGenerator
class SplineGenerator:
def __call__(self, INP):
us = self._un_sp
if not isinstance(INP, oofun):
raise FuncDesignerException('for scipy_InterpolatedUnivariateSpline input should be oovar/oofun,other cases not implemented yet')
def d(x):
x = np.asfarray(x)
#if x.size != 1:
#raise FuncDesignerException('for scipy_InterpolatedUnivariateSpline input should be oovar/oofun with output size = 1,other cases not implemented yet')
return us.__call__(x, 1)
def f(x):
x = np.asfarray(x)
#if x.size != 1:
#raise FuncDesignerException('for scipy_InterpolatedUnivariateSpline input should be oovar/oofun with output size = 1,other cases not implemented yet')
tmp = us.__call__(x.flatten() if x.ndim > 1 else x)
return tmp if x.ndim <= 1 else tmp.reshape(x.shape)
r = oofun(f, INP, d = d, isCostly=True, vectorized=True)
r._nonmonotone_x = self._nonmonotone_x
r._nonmonotone_y = self._nonmonotone_y
diffX, diffY = np.diff(self._X), np.diff(self._Y)
if (all(diffX >= 0) or all(diffX <= 0)) and (all(diffY >= 0) or all(diffY <= 0)) and self._k in (1, 3):
r.criticalPoints = False
elif self._k == 1:
r._interval = lambda *args: spline_interval_analysis_engine(r, *args)
else:
def _interval(*args, **kw):
raise FuncDesignerException('''
Currently interval calculations are implemented for
sorted monotone splines with order 1 or 3 only''')
r._interval = _interval
def Plot():
print('Warning! Plotting spline is recommended from FD spline generator, not initialized spline')
self.plot()
def Residual():
print('Warning! Getting spline residual is recommended from FD spline generator, not initialized spline')
return self.residual()
r.plot, r.residual = Plot, Residual
return r
def __init__(self, us, *args, **kwargs):
self._un_sp = us
self._X, self._Y = np.asfarray(args[0]), np.asfarray(args[1])
diffY = np.diff(self._Y)
ind_nonmonotone = np.where(diffY[1:] * diffY[:-1] < 0)[0] + 1
self._nonmonotone_x = self._X[ind_nonmonotone]
self._nonmonotone_y = self._Y[ind_nonmonotone]
if len(args) >= 5:
k = args[4]
elif 'k' in kwargs:
k = kwargs['k']
else:
k = 3 # default for InterpolatedUnivariateSpline
self._k = k
def plot(self):
try:
import pylab
except:
print('You should have matplotlib installed')
return
pylab.scatter(self._X, self._Y, marker='o')
YY = self._un_sp.__call__(self._X)
pylab.plot(self._X, YY)
pylab.grid('on')
pylab.title('FuncDesigner spline checker')
pylab.show()
def residual(self):
YY = self._un_sp.__call__(self._X)
return np.max(np.abs(YY - self._Y))
def spline_interval_analysis_engine(S, domain, dtype):
lb_ub, definiteRange = S.input[0]._interval(domain, dtype)
lb, ub = lb_ub[0], lb_ub[1]
x, y = S._nonmonotone_x, S._nonmonotone_y
tmp = S.fun(lb_ub)
_inf, _sup = tmp[0], tmp[1]
for i, xx in enumerate(x):
yy = y[i]
ind = np.logical_and(lb < xx, xx < ub)
_inf[ind] = np.where(_inf[ind] < yy, _inf[ind], yy)
_sup[ind] = np.where(_sup[ind] > yy, _sup[ind], yy)
r = np.vstack((_inf, _sup))
# TODO: modify definiteRange for out-of-bounds splines
# definiteRange = False
return r, definiteRange
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/Interval.py 000644 000766 000024 00000015705 12100431107 025553 0 ustar 00rfadams staff 000000 000000 from numpy import ndarray, asscalar, isscalar, floor, pi, \
copy as Copy, logical_and, where, asarray, any, atleast_1d, vstack, \
searchsorted
import numpy as np
try:
from bottleneck import nanmin, nanmax
except ImportError:
from numpy import nanmin, nanmax
class Interval:
def __init__(self, l, u, definiteRange):
if isinstance(l, ndarray) and l.size == 1: l = asscalar(l)
if isinstance(u, ndarray) and u.size == 1: u = asscalar(u)
self.lb, self.ub, self.definiteRange = l, u, definiteRange
def __str__(self):
return 'FuncDesigner interval with lower bound %s and upper bound %s' % (self.lb, self.ub)
def __repr__(self):
return str(self)
def ZeroCriticalPoints(lb_ub):
arg_infinum, arg_supremum = lb_ub[0], lb_ub[1]
if isscalar(arg_infinum):
return [0.0] if arg_infinum < 0.0 < arg_supremum else []
tmp = Copy(arg_infinum)
#tmp[where(logical_and(arg_infinum < 0.0, arg_supremum > 0.0))] = 0.0
tmp[atleast_1d(logical_and(arg_infinum < 0.0, arg_supremum > 0.0))] = 0.0
return [tmp]
#def IntegerCriticalPoints(arg_infinum, arg_supremum):
# # TODO: check it for rounding errors
# return arange(ceil(arg_infinum), ceil(1.0+arg_supremum), dtype=float).tolist()
#def TrigonometryCriticalPoints2(arg_infinum, arg_supremum):
# n1, n2 = int(floor(2 * arg_infinum / pi)), int(ceil(2 * arg_supremum / pi))
# # 6 instead of 5 for more safety, e.g. small numerical rounding effects
# return [i / 2.0 * pi for i in range(n1, amin((n1+6, n2))) if (arg_infinum < (i / 2.0) * pi < arg_supremum)]
# TODO: split TrigonometryCriticalPoints into (pi/2) *(2k+1) and (pi/2) *(2k)
def TrigonometryCriticalPoints(lb_ub):
arg_infinum, arg_supremum = lb_ub[0], lb_ub[1]
# returns points with coords n * pi/2, arg_infinum <= n * pi/2<= arg_supremum,n -array of integers
arrN = asarray(atleast_1d(floor(2 * arg_infinum / pi)), int)
Tmp = []
for i in range(1, 6):
th = (arrN+i)*pi/2
#ind = where(logical_and(arg_infinum < th, th < arg_supremum))[0]
ind = logical_and(arg_infinum < th, th < arg_supremum)
#if ind.size == 0: break
if not any(ind): break
tmp = atleast_1d(Copy(arg_infinum))
tmp[atleast_1d(ind)] = asarray((arrN[ind]+i)*pi/2, dtype = tmp.dtype)
Tmp.append(tmp)
return Tmp
# 6 instead of 5 for more safety, e.g. small numerical rounding effects
#return [i / 2.0 * pi for i in range(n1, amin((n1+6, n2))) if (arg_infinum < (i / 2.0) * pi < arg_supremum)]
#def halph_pi_x_2k_plus_one_points(arg_infinum, arg_supremum):
# n1 = asarray(floor(2 * arg_infinum / pi), int)
# Tmp = []
# for i in range(1, 7):
# if i% 2: continue
# ind = where(logical_and(arg_infinum < (n1+i)*pi/2, (n1+i)*pi/2< arg_supremum))[0]
# if ind.size == 0: break
# tmp = arg_infinum.copy()
# #assert (n1+i)*pi/2 < 6.3
# tmp[ind] = (n1[ind]+i)*pi/2
# Tmp.append(tmp)
# #raise 0
# return Tmp
#
def ZeroCriticalPointsInterval(inp, func):
def interval(domain, dtype):
lb_ub, definiteRange = inp._interval(domain, dtype)
lb, ub = lb_ub[0], lb_ub[1]
ind1, ind2 = lb < 0.0, ub > 0.0
ind = logical_and(ind1, ind2)
tmp = vstack((lb, ub))
TMP = func(tmp)
t_min, t_max = atleast_1d(nanmin(TMP, 0)), atleast_1d(nanmax(TMP, 0))
if any(ind):
F0 = func(0.0)
t_min[atleast_1d(logical_and(ind, t_min > F0))] = F0
t_max[atleast_1d(logical_and(ind, t_max < F0))] = F0
return vstack((t_min, t_max)), definiteRange
return interval
# if isscalar(arg_infinum):
# return [0.0] if arg_infinum < 0.0 < arg_supremum else []
# tmp = Copy(arg_infinum)
# #tmp[where(logical_and(arg_infinum < 0.0, arg_supremum > 0.0))] = 0.0
# tmp[atleast_1d(logical_and(arg_infinum < 0.0, arg_supremum > 0.0))] = 0.0
# return [tmp]
def nonnegative_interval(inp, func, domain, dtype, F0, shift = 0.0):
lb_ub, definiteRange = inp._interval(domain, dtype)
lb, ub = lb_ub[0], lb_ub[1]
t_min_max = func(lb_ub)
th = shift # 0.0 + shift = shift
ind = lb < th
if any(ind):
t_min_max[0][atleast_1d(logical_and(lb < th, ub >= th))] = F0
if definiteRange is not False:
if definiteRange is True:
definiteRange = np.empty_like(lb)
definiteRange.fill(True)
definiteRange[ind] = False
return t_min_max, definiteRange
def box_1_interval(inp, func, domain, dtype, F_l, F_u):
lb_ub, definiteRange = inp._interval(domain, dtype)
lb, ub = lb_ub[0], lb_ub[1]
t_min_max = func(lb_ub)
ind = lb < -1
if any(ind):
t_min_max[0][atleast_1d(logical_and(lb < -1, ub >= -1))] = F_l
if definiteRange is not False:
if definiteRange is True:
definiteRange = np.empty_like(lb)
definiteRange.fill(True)
definiteRange[ind] = False
ind = ub > 1
if any(ind):
t_min_max[0][atleast_1d(logical_and(lb < 1, ub >= 1))] = F_u
if definiteRange is not False:
if definiteRange is True:
definiteRange = np.empty_like(lb)
definiteRange.fill(True)
definiteRange[ind] = False
return t_min_max, definiteRange
# lb, ub = lb_ub[0], lb_ub[1]
# ind = lb < -1.0
# if any(ind):
# t_min, t_max = atleast_1d(empty_like(lb)), atleast_1d(empty_like(ub))
# t_min.fill(nan)
# t_max.fill(nan)
# ind2 = ub >= -1.0
# t_min[atleast_1d(logical_not(ind))] = func(lb)
# t_min[atleast_1d(logical_and(ind2, ind))] = 0.0
# t_max[atleast_1d(ind2)] = func(ub[ind2])
# t_min, t_max = func(lb), func(ub)
# t_min[atleast_1d(logical_and(lb < 0, ub > 0))] = 0.0
#
# # TODO: rework it with matrix operations
# definiteRange = False
#
# else:
# t_min, t_max = func(lb), func(ub)
# return vstack((t_min, t_max)), definiteRange
def adjust_lx_WithDiscreteDomain(Lx, v):
if v.domain is bool or v.domain is 'bool':
Lx[Lx != 0] = 1
else:
d = v.domain
ind = searchsorted(d, Lx, 'left')
ind2 = searchsorted(d, Lx, 'right')
ind3 = where(ind!=ind2)[0]
#Tmp = Lx[:, ind3].copy()
Tmp = d[ind[ind3]]
#if any(ind==d.size):print 'asdf'
ind[ind==d.size] -= 1# Is it ever encountered?
# ind[ind==d.size-1] -= 1
Lx[:] = d[ind]
Lx[ind3] = asarray(Tmp, dtype=Lx.dtype)
def adjust_ux_WithDiscreteDomain(Ux, v):
if v.domain is bool or v.domain is 'bool':
Ux[Ux != 1] = 0
else:
d = v.domain
ind = searchsorted(d, Ux, 'left')
ind2 = searchsorted(d, Ux, 'right')
ind3 = where(ind!=ind2)[0]
#Tmp = Ux[:, ind3].copy()
Tmp = d[ind[ind3]]
#ind[ind==d.size] -= 1
ind[ind==0] = 1
Ux[:] = d[ind-1]
Ux[ind3] = asarray(Tmp, dtype=Ux.dtype)
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/logic.py 000644 000766 000024 00000000242 12100431107 025052 0 ustar 00rfadams staff 000000 000000 from FDmisc import FuncDesignerException
from .ooFun import oofun, BooleanOOFun, nlh_and, nlh_not
import numpy as np
#def NOT(Args):
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/multiarray.py 000644 000766 000024 00000006643 12100431107 026161 0 ustar 00rfadams staff 000000 000000 import numpy as np, operator
from numpy import ndarray
from FuncDesigner.FDmisc import FuncDesignerException
from baseClasses import MultiArray
if 'div' in operator.__dict__:
div = operator.div
else:
div = operator.truediv
# TODO: rework buggy multiarray.size
#delattr(multiarray, 'size')
class multiarray(MultiArray):
__array_priority__ = 5
__add__ = lambda self, other: multiarray_op(self, other, operator.add)
__radd__ = lambda self, other: self.__add__(other)
__sub__ = lambda self, other: multiarray_op(self, other, operator.sub)
__rsub__ = lambda self, other: multiarray_op(-self, other, operator.add)
__mul__ = lambda self, other: multiarray_op(self, other, operator.mul)
__rmul__ = lambda self, other: self.__mul__(other)
__div__ = lambda self, other: multiarray_op(self, other, div)
__rdiv__ = lambda self, other: multiarray_op(other, self, div)
__truediv__ = __div__
__rtruediv__ = __rdiv__
# TODO: rdiv, rpow
__pow__ = lambda self, other: multiarray_op(self, other, operator.pow)
__rpow__ = lambda self, other: multiarray_op(other, self, operator.pow)
def __getitem__(self, ind):
# if self.ndim <= 1:
# if ind is 0:
# return self
# else:
# print('multiarray ind:', ind)
return self.view(np.ndarray)[:, ind].view(multiarray) if type(ind) in (int, np.int32, np.int64, np.int16, np.int8) \
else self.__getslice__(ind.start, ind.stop) if type(ind) != tuple \
else self.__getslice__(ind[0], ind[1])
def __getslice__(self, ind1, ind2):
#TODO: mb check if size is known then use it instead of None?
cond_1 = ind1 is None
cond_2 = ind2 is None or (type(ind2) == slice and ind2.start is None and ind2.stop is None and ind2.step is None)
if cond_1 and cond_2:
return self
if cond_1:
ind1 = 0
if cond_2:
ind2 = self.shape[1]
return self.view(np.ndarray)[:, ind1:ind2].view(multiarray)
# TODO: check it!
#toarray = lambda self: self.view(ndarray)
def sum(self, *args, **kw):
if any([v is not None for v in args]): # somehow triggered from pswarm
raise FuncDesignerException('arguments for FD multiarray sum are not implemened yet')
if any([v is not None for v in kw.values()]):
raise FuncDesignerException('keyword arguments for FD multiarray sum are not implemened yet')
tmp = self.reshape(-1, 1) if self.ndim < 2 else self
return np.sum(tmp.view(ndarray), 1).view(multiarray)
def multiarray_op(x, y, op):
if isinstance(y, multiarray):
Y = y.reshape(-1, 1) if y.ndim < 2 else y
if isinstance(x, multiarray):
assert x.ndim < 3 and y.ndim < 3, 'unimplemented yet'
X = x.reshape(-1, 1) if x.ndim < 2 else x
r = op(X.view(ndarray), Y.view(ndarray))
else:
r = op(x.reshape(-1, 1) if isinstance(x, ndarray) and x.size != 1 else x, Y.view(ndarray).T).T
elif isinstance(x, multiarray): # and y is not multiarray here
X = x.reshape(-1, 1) if x.ndim < 2 else x
r = op(X.view(ndarray), y.reshape(1, -1) if isinstance(y, ndarray) and y.size != 1 else y)
else: # neither x nor y are multiarrays
raise FuncDesignerException('bug in FuncDesigner kernel')
return r.view(multiarray)
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/ode.py 000644 000766 000024 00000020554 12100431107 024534 0 ustar 00rfadams staff 000000 000000 from translator import FuncDesignerTranslator
from FDmisc import FuncDesignerException, _getDiffVarsID
from numpy import ndarray, hstack, vstack, isscalar, asarray
from ooVar import oovar
from ooFun import atleast_oofun
class ode:
# Ordinary differential equations
_isInitialized = False
solver = 'scipy_lsoda'
def __init__(self, equations, startPoint, *args, **kwargs):
if len(args) > 2:
raise FuncDesignerException('incorrect ode definition, too many args are obtained')
if not isinstance(equations, dict):
raise FuncDesignerException('1st argument of ode constructor should be Python dict')
if not isinstance(startPoint, dict):
raise FuncDesignerException('2nd argument of ode constructor should be Python dict')
if len(args) == 2:
print("""
FuncDesigner warning: you're using obsolete ode() API,
now time should be passed as 3rd argument {timeVariable: times}
or just array of times if equations are time-independend""")
timeVariable, times = args[0], args[1]
else:
if len(args) != 1 or type(args[0]) not in (dict, list, tuple, ndarray):
raise FuncDesignerException('3rd argument of ode constructor must be dict {timeVariable:array_of_times} or just array_of_times')
if type(args[0]) == dict:
if len(args[0]) != 1:
raise FuncDesignerException('in time argument dict has to have exactly 1 entry')
timeVariable, times = list(args[0].items())[0]
if not isinstance(timeVariable, oovar):
raise FuncDesignerException('incorrect time variable, must be FuncDesigner oofun')
else:
times = args[0]
timeVariable = None
self.timeVariable = timeVariable
if timeVariable is not None and timeVariable in equations:
raise FuncDesignerException("ode: differentiation of a variable by itself (time by time) is treated as a potential bug and thus is forbidden")
if not (isinstance(times, (list, tuple)) or (isinstance(times, ndarray) and times.ndim == 1)):
raise FuncDesignerException('incorrect user-defined time argument, must be Python list or numpy array of time values')
self.times = times
self._fd_func, self._startPoint, self._timeVariable, self._times, self._kwargs = equations, startPoint, timeVariable, times, kwargs
startPoint = dict([(key, val) for key, val in startPoint.items()])
if timeVariable is not None:
startPoint[timeVariable] = times[0]
y0 = []
Funcs = []
# setting oovar.size is risky - it can affect code written after the ode is solved
# thus Point4TranslatorAssignment is used instead
Point4TranslatorAssignment = {}
for v, func in equations.items():
func = atleast_oofun(func)
# if not isinstance(func, oofun):
# func = fixed_oofun(func)
if not isinstance(v, oovar):
raise FuncDesignerException('ode: dict keys must be FuncDesigner oovars, got "%s" instead' % type(v))
startFVal = asarray(func(startPoint))
y0.append(asarray(startPoint[v]))
Funcs.append(func)
Point4TranslatorAssignment[v] = startFVal
if startFVal.size != asarray(startPoint[v]).size or (hasattr(v, 'size') and isscalar(v.size) and startFVal.size != v.size):
raise FuncDesignerException('error in user-defined data: oovar "%s" size is not equal to related function value in start point' % v.name)
self.y0 = hstack(y0)
self.varSizes = [y.size for y in y0]
ooT = FuncDesignerTranslator(Point4TranslatorAssignment)
self.ooT = ooT
def func (y, t):
tmp = dict(ooT.vector2point(y))
if timeVariable is not None:
tmp[timeVariable] = t
return hstack([func(tmp) for func in Funcs])
self.func = func
_FDVarsID = _getDiffVarsID()
def derivative(y, t):
tmp = dict(ooT.vector2point(y))
if timeVariable is not None:
tmp[timeVariable] = t
r = []
for func in Funcs:
tt = func.D(tmp, fixedVarsScheduleID = _FDVarsID)
if timeVariable is not None:
tt.pop(timeVariable)
r.append(ooT.pointDerivative2array(tt))
return vstack(r)
self.derivative = derivative
self.Point4TranslatorAssignment = Point4TranslatorAssignment
def solve(self, solver='scipy_lsoda', *args, **kwargs): # mb for future implementation - add some **kwargs here as well
if len(args) > 0:
raise FuncDesignerException('no args are currently available for the function ode::solve')
solverName = solver if isinstance(solver, str) else solver.__name__
if solverName.startswith('interalg'):
try:
from openopt import ODE
except ImportError:
raise FuncDesignerException('You should have openopt insalled')
prob = ODE(self._fd_func, self._startPoint, **self._kwargs)
prob.timeVariable = self.timeVariable
prob.times = self.times
r = prob.solve(solver, **kwargs)
y_var = list(prob._x0.keys())[0]
res = 0.5 * (prob.extras[y_var]['infinums'] + prob.extras[y_var]['supremums'])
times = 0.5 * (prob.extras['startTimes'] + prob.extras['endTimes'])
if len(self._times) != 2:
# old
from scipy.interpolate import UnivariateSpline
if 'ftol' in kwargs.keys():
s = self._kwargs['ftol']
elif 'fTol' in kwargs.keys():
s = self._kwargs['fTol']
elif 'ftol' in self._kwargs.keys():
s = self._kwargs['ftol']
elif 'fTol' in self._kwargs.keys():
s = self._kwargs['fTol']
# walkaround a bug with essential slowdown
if times[-1] < times[0]:
times = times[::-1]
res = res[::-1]
interp = UnivariateSpline(times, res, k=1, s=s**2)
times = self._times
res = interp(times)
# new
# ind = searchsorted(prob.extras['startTimes'], self._times, 'right')
# ind[ind==times.size] = times.size-1
# tmp = res[ind]
## lti, rti = self._times - prob.extras['startTimes'][ind], prob.extras['endTimes'][ind] - self._times
## tmp = (lti * res[ind] + rti * res[ind]) / (lti + rti)
## tmp[logical_and(lti==0, rti==0)] = res[ind]
# res = tmp
r.xf = {y_var:res}
r._xf = {y_var.name: res}
if self.timeVariable is not None:
r.xf[self._timeVariable] = times
r._xf[self._timeVariable.name] = times
else:
if solver != 'scipy_lsoda': raise FuncDesignerException('incorrect ODE solver')
try:
from scipy import integrate
except ImportError:
raise FuncDesignerException('to solve ode you mush have scipy installed, see http://openop.org/SciPy')
y, infodict = integrate.odeint(self.func, self.y0, self.times, Dfun = self.derivative, full_output=True)
resultDict = dict(self.ooT.vector2point(y.T))
for key, value in resultDict.items():
if min(value.shape) == 1:
resultDict[key] = value.flatten()
r = FuncDesigner_ODE_Result(resultDict)
r.msg = infodict['message']
r.extras = {'infodict': infodict}
return r
class FuncDesigner_ODE_Result:
# TODO: prevent code clone with runprobsolver.py
def __init__(self, resultDict):
self.xf = resultDict
if not hasattr(self, '_xf'):
self._xf = dict([(var.name, value) for var, value in resultDict.items()])
def __call__(self, *args):
r = [(self._xf[arg] if isinstance(arg, str) else self.xf[arg]) for arg in args]
return r[0] if len(args)==1 else r
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/ooarray.py 000644 000766 000024 00000032663 12100431107 025445 0 ustar 00rfadams staff 000000 000000 from baseClasses import OOArray
from FuncDesigner.multiarray import multiarray
from ooFun import oofun, Constraint
from numpy import isscalar, asscalar, ndarray, asarray, atleast_1d, asanyarray
import numpy as np
from FDmisc import FuncDesignerException
class ooarray(OOArray):
__array_priority__ = 25 # !!! it should exceed oofun.__array_priority__ !!!
_is_array_of_oovars = False
def __new__(self, *args, **kwargs):
#assert len(kwargs) == 0
tmp = args[0] if len(args) == 1 else args
obj = asarray(tmp).view(self)
#if obj.ndim != 1: raise FuncDesignerException('only 1-d ooarrays are implemented now')
#if obj.dtype != object:obj = np.asfarray(obj) #TODO: FIXME !
obj._id = oofun._id
obj.name = 'unnamed_ooarray_%d' % obj._id
oofun._id += 1
return obj
# def __init__(self, *args, **kw):
# self._id = oofun._id
# self.name = 'unnamed_ooarray_%d' % self._id
# oofun._id += 1
__hash__ = lambda self: self._id
def __len__(self):
return self.size
expected_kwargs = set(('tol', 'name'))
def __call__(self, *args, **kwargs):
#if self.dtype != object: return self.view(ndarray)
# TODO: give different names for each element while assigning name to ooarray
expected_kwargs = self.expected_kwargs
#if not set(kwargs.keys()).issubset(expected_kwargs):
#raise FuncDesignerException('Unexpected kwargs: should be in '+str(expected_kwargs)+' got: '+str(kwargs.keys()))
for elem in expected_kwargs:
if elem in kwargs:
setattr(self, elem, kwargs[elem])
if len(args) > 1: raise FuncDesignerException('No more than single argument is expected')
if len(args) == 0:
if len(kwargs) == 0: raise FuncDesignerException('You should provide at least one argument')
#return self
if len(args) != 0 and isinstance(args[0], str):
self.name = args[0]
for i, elem in enumerate(self.view(ndarray)):
if isinstance(elem, oofun):
elem(self.name + '_' + str(i))
args = args[1:]
if len(args) == 0:
return self
#tmp = asarray([asscalar(asarray(self[i](*args, **kwargs))) if isinstance(self[i], oofun) else self[i] for i in range(self.size)])
# TODO: get rid of self in args[0]
if self._is_array_of_oovars and isinstance(args[0], dict) and self in args[0] and len(args) == 1 and len(kwargs) == 0:
return args[0][self]
Tmp = [self[i](*args, **kwargs) if isinstance(self[i], oofun) else self[i] for i in range(self.size)]
tmp = asanyarray(Tmp)
if np.any([isinstance(elem, multiarray) for elem in Tmp]):
tmp = tmp.T.view(multiarray)
if tmp.ndim == 2 or tmp.dtype != object:
return tmp
else:
#tmp = tmp.flatten()
return ooarray(tmp)
def __mul__(self, other):
if self.size == 1:
return ooarray(asscalar(self)*other)
elif isscalar(other):
# TODO: mb return mere ooarray(self.view(ndarray)*other) or other.view(ndarray)
return ooarray(self.view(ndarray)*other if self.dtype != object else [self[i]*other for i in range(self.size)])
elif isinstance(other, oofun):
hasSize = 'size' in dir(other)
if not hasSize:
# print('''
# FuncDesigner warning:
# to perform the operation
# (ooarray multiplication on oofun)
# oofun size should be known.
# Assuming oofun size is 1,
# the value is ascribed to the oofun attributes.
# Handling of the issue is intended to be
# enhanced in future.''')
other.size = 1
#raise FuncDesignerException('to perform the operation oofun size should be known')
if other.size == 1:
if any([isinstance(elem, oofun) for elem in atleast_1d(self)]):
#if self.dtype == object:
s = atleast_1d(self)
return ooarray([s[i]*other for i in range(self.size)])
else:
return ooarray(self*other)
else: # other.size > 1
# and self.size != 1
s, o = atleast_1d(self), atleast_1d(other)
return ooarray([s[i]*o[i] for i in range(self.size)])
elif isinstance(other, ndarray):
# TODO: mb return mere ooarray(self.view(ndarray)*other)?or other.view(ndarray)
return ooarray(self*asscalar(other) if other.size == 1 else [self[i]*other[i] for i in range(other.size)])
elif type(other) in (list, tuple):
r = self * asarray(other)
return r
else:
raise FuncDesignerException('bug in multiplication')
def __div__(self, other):
if self.size == 1:
return asscalar(self)/other
elif isscalar(other) or (isinstance(other, ndarray) and other.size in (1, self.size)):
return self * (1.0/other)
elif isinstance(other, oofun):
if self.dtype != object:
return self.view(ndarray) / other
else:
s = atleast_1d(self)
return ooarray([s[i] / other for i in range(self.size)])
elif isinstance(other, ooarray):
if self.dtype != object:
return self.view(ndarray) / other.view(ndarray)
else:
# TODO: mb return mere ooarray(self.view(ndarray) / other)? or other.view(ndarray)
s, o = atleast_1d(self), atleast_1d(other)
return ooarray([s[i] / o[i] for i in range(self.size)])
else:
raise FuncDesignerException('unimplemented yet')
__truediv__ = __div__
__floordiv__ = __div__
def __rdiv__(self, other):
if self.size == 1:
return other / asscalar(self)
return ooarray([1.0 / elem for elem in self.view(ndarray)]) * other
__rtruediv__ = __rdiv__
def __add__(self, other):
if isinstance(other, list):
other = ooarray(other)
if isscalar(other) or (isinstance(other, ndarray) and other.size in (1, self.size)):
r = ooarray(self.view(ndarray) + other)
elif isinstance(other, oofun):
if self.dtype != object:
r = self.view(ndarray) + other
else:
s = atleast_1d(self)
r = ooarray([s[i] + other for i in range(self.size)])
elif isinstance(other, ndarray):
if self.dtype != object:
r = self.view(ndarray) + other.view(ndarray)
elif self.size == 1:
r = other + asscalar(self)
else:
# TODO: mb return mere ooarray(self.view(ndarray) + other) or ooarray(self.view(ndarray) + other.view(ndarray))?
r = ooarray([self[i] + other[i] for i in range(self.size)])
else:
raise FuncDesignerException('unimplemented yet')
if isinstance(r, ndarray) and r.size == 1:
r = asscalar(r)
return r
# TODO: check why it doesn't work with oofuns
def __radd__(self, other):
return self + other
def __rmul__(self, other):
return self * other
# TODO : fix it
# def __rdiv__(self, other):
# return self * other
def __pow__(self, other):
if isinstance(other, ndarray) and other.size > 1 and self.size > 1:
return ooarray([self[i]**other[i] for i in range(self.size)])
Self = atleast_1d(self.view(ndarray))
if any([isinstance(elem, oofun) for elem in Self]):
#if self.dtype == object:
return ooarray([elem**other for elem in Self])
# TODO: is this part of code trigger any time?
return self.view(ndarray)**other
def __rpow__(self, other):
if isscalar(other) or ('size' in dir(other) and isscalar(other.size) and other.size == 1) or 'size' not in dir(other):
return ooarray([other ** elem for elem in self.tolist()])
return ooarray([other[i] ** elem for i, elem in enumerate(self.tolist())])
def __eq__(self, other):
r = self - other
if r.dtype != object: return all(r)
if r.size == 1: return asscalar(r)==0
# TODO: rework it
return ooarray([Constraint(elem, lb=0.0, ub=0.0) for elem in r.tolist()])
#else: raise FuncDesignerException('unimplemented yet')
def __lt__(self, other):
if self.dtype != object and (not isinstance(other, ooarray) or other.dtype != object):
return self.view(ndarray) < (other.view(ndarray) if isinstance(other, ooarray) else other)
if isinstance(other, (ndarray, list, tuple)) and self.size > 1 and len(other) > 1:
return ooarray([self[i] < other[i] for i in range(self.size)])
if isscalar(other) or (isinstance(other, (ndarray, list, tuple)) and len(other) == 1):
return ooarray([elem < other for elem in self])
if isinstance(other, oofun):
if 'size' in other.__dict__ and not isinstance(other.size, oofun):
if other.size == self.size:
return ooarray([elem[i] < other[i] for i in range(self.size)])
elif self.size == 1:
return ooarray([self[0] < other[i] for i in range(other.size)])
else:
FuncDesignerException('bug or yet unimplemented case in FD kernel')
else:
# !!! assunimg other.size = 1
return ooarray([elem < other for elem in self])
raise FuncDesignerException('unimplemented yet')
def __le__(self, other):
if self.dtype != object and (not isinstance(other, ooarray) or other.dtype != object):
return self.view(ndarray) <= (other.view(ndarray) if isinstance(other, ooarray) else other)
if isinstance(other, (ndarray, list, tuple)) and self.size > 1 and len(other) > 1:
return ooarray([self[i] <= other[i] for i in range(self.size)])
if isscalar(other) or (isinstance(other, (ndarray, list, tuple)) and len(other) == 1):
return ooarray([elem <= other for elem in self])
if isinstance(other, oofun):
if 'size' in other.__dict__ and not isinstance(other.size, oofun):
if other.size == self.size:
return ooarray([elem[i] <= other[i] for i in range(self.size)])
elif self.size == 1:
return ooarray([self[0] <= other[i] for i in range(other.size)])
else:
FuncDesignerException('bug or yet unimplemented case in FD kernel')
else:
# !!! assunimg other.size = 1
return ooarray([elem <= other for elem in self])
raise FuncDesignerException('unimplemented yet')
def __gt__(self, other):
if self.dtype != object and (not isinstance(other, ooarray) or other.dtype != object):
return self.view(ndarray) > (other.view(ndarray) if isinstance(other, ooarray) else other)
if isinstance(other, (ndarray, list, tuple)) and self.size > 1 and len(other) > 1:
return ooarray([self[i] > other[i] for i in range(self.size)])
if isscalar(other) or (isinstance(other, (ndarray, list, tuple)) and len(other) == 1):
return ooarray([elem > other for elem in self])
if isinstance(other, oofun):
if 'size' in other.__dict__ and not isinstance(other.size, oofun):
if other.size == self.size:
return ooarray([elem[i] > other[i] for i in range(self.size)])
elif self.size == 1:
return ooarray([self[0] > other[i] for i in range(other.size)])
else:
FuncDesignerException('bug or yet unimplemented case in FD kernel')
else:
# !!! assunimg other.size = 1
return ooarray([elem > other for elem in self])
raise FuncDesignerException('unimplemented yet')
def __ge__(self, other):
if self.dtype != object and (not isinstance(other, ooarray) or other.dtype != object):
return self.view(ndarray) >= (other.view(ndarray) if isinstance(other, ooarray) else other)
if isinstance(other, (ndarray, list, tuple)) and self.size > 1 and len(other) > 1:
return ooarray([self[i] >= other[i] for i in range(self.size)])
if isscalar(other) or (isinstance(other, (ndarray, list, tuple)) and len(other) == 1):
return ooarray([elem >= other for elem in self])
if isinstance(other, oofun):
if 'size' in other.__dict__ and not isinstance(other.size, oofun):
if other.size == self.size:
return ooarray([elem[i] >= other[i] for i in range(self.size)])
elif self.size == 1:
return ooarray([self[0] >= other[i] for i in range(other.size)])
else:
FuncDesignerException('bug or yet unimplemented case in FD kernel')
else:
# !!! assunimg other.size = 1
return ooarray([elem >= other for elem in self])
raise FuncDesignerException('unimplemented yet')
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/ooFun.py 000644 000766 000024 00000340677 12100431107 025066 0 ustar 00rfadams staff 000000 000000 # created by Dmitrey
PythonSum = sum
from numpy import inf, asfarray, copy, all, any, atleast_2d, zeros, dot, asarray, atleast_1d, \
ones, ndarray, where, array, nan, vstack, eye, array_equal, isscalar, log, hstack, sum as npSum, prod, nonzero,\
isnan, asscalar, zeros_like, ones_like, amin, amax, logical_and, logical_or, isinf, logical_not, logical_xor, flipud, \
tile, float64, searchsorted, int8, int16, int32, int64, isfinite, log2, string_, asanyarray, bool_
#from logic import AND
#from traceback import extract_stack
try:
from bottleneck import nanmin, nanmax
except ImportError:
from numpy import nanmin, nanmax
from FDmisc import FuncDesignerException, Diag, Eye, pWarn, scipyAbsentMsg, scipyInstalled, \
raise_except, DiagonalType, isPyPy
from ooPoint import ooPoint
from FuncDesigner.multiarray import multiarray
from Interval import Interval, adjust_lx_WithDiscreteDomain, adjust_ux_WithDiscreteDomain
import inspect
from baseClasses import OOArray, Stochastic
Copy = lambda arg: asscalar(arg) if type(arg)==ndarray and arg.size == 1 else arg.copy() if hasattr(arg, 'copy') else copy(arg)
Len = lambda x: 1 if isscalar(x) else x.size if type(x)==ndarray else x.values.size if isinstance(x, Stochastic) else len(x)
try:
from DerApproximator import get_d1, check_d1
DerApproximatorIsInstalled = True
except:
DerApproximatorIsInstalled = False
try:
import scipy
#from scipy import sparse
from scipy.sparse import hstack as HstackSP, vstack as VstackSP, isspmatrix_csc, isspmatrix_csr, eye as SP_eye, lil_matrix as SparseMatrixConstructor
def Hstack(Tuple):
ind = where([isscalar(elem) or prod(elem.shape)!=0 for elem in Tuple])[0].tolist()
elems = [Tuple[i] for i in ind]
if any([isspmatrix(elem) for elem in elems]):
return HstackSP(elems)
s = set([(0 if isscalar(elem) else elem.ndim) for elem in elems])
ndim = max(s)
if ndim <= 1: return hstack(elems)
assert ndim <= 2 and 1 not in s, 'bug in FuncDesigner kernel, inform developers'
return hstack(elems) if 0 not in s else hstack([atleast_2d(elem) for elem in elems])
def Vstack(Tuple):
ind = where([isscalar(elem) or prod(elem.shape)!=0 for elem in Tuple])[0].tolist()
elems = [Tuple[i] for i in ind]
if any([isspmatrix(elem) for elem in elems]):
return VstackSP(elems)
else:
return vstack(elems)
# s = set([(0 if isscalar(elem) else elem.ndim) for elem in elems])
# ndim = max(s)
# if ndim <= 1: return hstack(elems)
# assert ndim <= 2 and 1 not in s, 'bug in FuncDesigner kernel, inform developers'
# return hstack(elems) if 0 not in s else hstack([atleast_2d(elem) for elem in elems])
from scipy.sparse import isspmatrix
except:
scipy = None
isspmatrix = lambda *args, **kwargs: False
Hstack = hstack
Vstack = vstack
class oofun:
#TODO:
#is_oovarSlice = False
tol = 0.0
d = None # derivative
input = None#[]
#usedIn = set()
is_oovar = False
#is_stoch = False
isConstraint = False
#isDifferentiable = True
discrete = False
_isSum = False
_isProd = False
stencil = 3 # used for DerApproximator
#TODO: modify for cases where output can be partial
evals = 0
same = 0
same_d = 0
evals_d = 0
# finite-difference aproximation step
diffInt = 1.5e-8
maxViolation = 1e-2
_unnamedFunNumber = 1
_lastDiffVarsID = 0
_lastFuncVarsID = 0
_lastOrderVarsID = 0
criticalPoints = None
vectorized = False
getDefiniteRange = None
_neg_elem = None # used in render into quadratic
_usedIn = 0
_level = 0
#_directlyDwasInwolved = False
_id = 0
_BroadCastID = 0
_broadcast_id = 0
_point_id = 0
_point_id1 = 0
_f_key_prev = None
_f_val_prev = None
_d_key_prev = None
_d_val_prev = None
__array_priority__ = 15# set it greater than 1 to prevent invoking numpy array __mul__ etc
hasDefiniteRange = True
pWarn = lambda self, msg: pWarn(msg)
def disp(self, msg):
print(msg)
nlh = lambda self, *args, **kw: raise_except("probably you have involved boolean operation on continuous function, that is error")
lh = lambda self, *args, **kw: raise_except("probably you have involved boolean operation on continuous function, that is error")
def __getattr__(self, attr):
if attr == '__len__':
# TODO: fix it
if isPyPy:
return 1
else:
raise AttributeError('using len(oofun) is not possible yet, try using oofun.size instead')
elif attr == 'isUncycled':
self._getDep()
return self.isUncycled
elif attr == 'isCostly':
return self.d is None
elif attr != 'size':
raise AttributeError('you are trying to obtain incorrect attribute "%s" for FuncDesigner oofun "%s"' %(attr, self.name))
# to prevent creating of several oofuns binded to same oofun.size
r = oofun(lambda x: asarray(x).size, self, discrete = True, getOrder = lambda *args, **kwargs: 0)
self.size = r
return r
""" Class constructor """
def __init__(self, fun, input=None, *args, **kwargs):
#def __init__(self, fun, input=[], *args, **kwargs):
assert len(args) == 0 #and input is not None
self.fun = fun
#self._broadcast_id = 0
self._id = oofun._id
self.attachedConstraints = set()
self.args = ()
oofun._id += 1 # CHECK: it should be int32! Other types cannot be has keys!
if 'name' not in kwargs.keys():
self.name = 'unnamed_oofun_' + str(oofun._unnamedFunNumber)
oofun._unnamedFunNumber += 1
for key, item in kwargs.items():
#assert key in self.__allowedFields__ # TODO: make set comparison
setattr(self, key, item)
if isinstance(input, (tuple, list)):
self.input = [(elem if isinstance(elem, (oofun, OOArray)) else array(elem, 'float')) for elem in input]
elif input is not None:
self.input = [input]
else:
self.input = [None] # TODO: get rid of None, use input = [] instead
# TODO: fix it for ooarray!
# if input is not None:
# levels = [0]
# for elem in self.input: # if a
# if isinstance(elem, oofun):
# elem._usedIn += 1
# levels.append(elem._level)
# self._level = max(levels)+1
__hash__ = lambda self: self._id
def attach(self, *args, **kwargs):
if len(kwargs) != 0:
raise FuncDesignerException('keyword arguments are not implemented for FuncDesigner function "attach"')
assert len(args) != 0
Args = args[0] if len(args) == 1 and type(args[0]) in (tuple, list, set) else args
for arg in Args:
if not isinstance(arg, BaseFDConstraint):
raise FuncDesignerException('the FD function "attach" currently expects only constraints')
self.attachedConstraints.update(Args)
return self
def removeAttachedConstraints(self):
self.attachedConstraints = set()
__repr__ = lambda self: self.name
def _interval_(self, domain, dtype):
criticalPointsFunc = self.criticalPoints
if len(self.input) == 1 and criticalPointsFunc is not None:
arg_lb_ub, definiteRange = self.input[0]._interval(domain, dtype)
arg_infinum, arg_supremum = arg_lb_ub[0], arg_lb_ub[1]
if (not isscalar(arg_infinum) and arg_infinum.size > 1) and not self.vectorized:
raise FuncDesignerException('not implemented for vectorized oovars yet')
tmp = [arg_lb_ub]
if criticalPointsFunc is not False:
tmp += criticalPointsFunc(arg_lb_ub)
Tmp = self.fun(vstack(tmp))
if self.getDefiniteRange is not None:
definiteRange = logical_and(definiteRange, self.getDefiniteRange(arg_infinum, arg_supremum))
if not self.hasDefiniteRange:
# TODO: rework it as matrix operations
definiteRange = False #logical_and(definiteRange, a)
#indDefinite = all(isfinite)
#Tmp = self.fun(array(tmp))
return vstack((nanmin(Tmp, 0), nanmax(Tmp, 0))), definiteRange
else:
raise FuncDesignerException('interval calculations are unimplemented for the oofun (%s) yet' % self.name)
def interval(self, domain, dtype = float, resetStoredIntervals = True):
if type(domain) != ooPoint:
domain = ooPoint(domain, skipArrayCast = True)
lb_ub, definiteRange = self._interval(domain, dtype) #if type(domain) != ooPoint else self._interval2(domain, dtype)
# TODO: MB GET RID OF IT?
if resetStoredIntervals:
domain.storedIntervals = {}
return Interval(lb_ub[0], lb_ub[1], definiteRange)
def _interval(self, domain, dtype):
tmp = domain.dictOfFixedFuncs.get(self, None)
if tmp is not None:
return tile(tmp, (2, 1)), True
v = domain.modificationVar
if v is None or ((v not in self._getDep() or self.is_oovar) and self is not v):
r = domain.storedIntervals.get(self, None)
if r is not None:
return r
if v is not None:
r = domain.localStoredIntervals.get(self, None)
if r is not None:
return r
r = self._interval_(domain, dtype)
if domain.useSave:
domain.storedIntervals[self] = r
if v is not None:
domain.localStoredIntervals[self] = r
return r
def iqg(self, domain, dtype = float, lb=None, ub=None, UB = None):
if type(domain) != ooPoint:
domain = ooPoint(domain, skipArrayCast=True)
domain.isMultiPoint=True
domain.useSave = True
r0 = self.interval(domain, dtype, resetStoredIntervals = False)
r0.lb, r0.ub = atleast_1d(r0.lb).copy(), atleast_1d(r0.ub).copy() # is copy required?
# TODO: get rid of useSave
domain.useSave = False
# TODO: rework it with indexation of required data
if lb is not None and ub is not None:
#debug start
# ind = logical_or(logical_or(r0.ub < lb, r0.lb > ub), all(logical_and(r0.lb >= lb, r0.ub <= ub)))
# sz = where(ind)[0]
# if sz.size != asarray(r0.lb).size:
# print('!', sz.size, asarray(r0.lb).size)
#debug end
ind = logical_or(logical_or(r0.ub < lb, r0.lb > ub), all(logical_and(r0.lb >= lb, r0.ub <= ub)))
elif UB is not None:
#print(where(r0.lb > UB)[0].size, asarray(r0.lb).size)
ind = r0.lb > UB
else:
ind = None
useSlicing = False
# print('-'*20)
# print(r0.lb.size)
if ind is not None:
if all(ind):
return {}, r0
j = where(~ind)[0]
#DOESN'T WORK FOR FIXED OOVARS AND DefiniteRange != TRUE YET
if 0 and j.size < 0.85*ind.size: # at least 15% of values to skip
useSlicing = True
tmp = []
for key, val in domain.storedIntervals.items():
Interval, definiteRange = val
if type(definiteRange) not in (bool, bool_):
definiteRange = definiteRange[j]
tmp.append((key, (Interval[:, j], definiteRange)))
_storedIntervals = dict(tmp)
Tmp = []
for key, val in domain.storedSums.items():
# TODO: rework it
R0, DefiniteRange0 = val.pop(-1)
#R0, DefiniteRange0 = val[-1]
R0 = R0[:, j]
if type(DefiniteRange0) not in (bool, bool_):
DefiniteRange0 = DefiniteRange0[j]
tmp = []
for k,v in val.items():
# TODO: rework it
# if k is (-1): continue
v = v[:, j]
tmp.append((k,v))
val = dict(tmp)
val[-1] = (R0, DefiniteRange0)
Tmp.append((key,val))
_storedSums = dict(Tmp)
#domain.storedSums = dict(tmp)
Tmp = []
for key, val in domain.items():
lb_,ub_ = val
# TODO: rework it when lb, ub will be implemented as 2-dimensional
Tmp.append((key, (lb_[j],ub_[j])))
dictOfFixedFuncs = domain.dictOfFixedFuncs
domain2 = ooPoint(Tmp, skipArrayCast=True)
domain2.storedSums = _storedSums
domain2.storedIntervals = _storedIntervals
domain2.dictOfFixedFuncs = dictOfFixedFuncs
domain2.isMultiPoint=True
domain = domain2
domain.useAsMutable = True
r = {}
Dep = (self._getDep() if not self.is_oovar else set([self])).intersection(domain.keys())
## #debug
# if useSlicing:
## print('ind', ind)
## print('r0.ub:', r0.ub)
# r1 = {}
# for i, v in enumerate(Dep):
# domain.modificationVar = v
# r_l, r_u = self._iqg(domain, dtype, r0)
#
# domain2.useAsMutable = True
# domain2.modificationVar = v
# r_l2, r_u2 = self._iqg(domain2, dtype, r0)
# if useSlicing and not (r_l is r0):# r_l is r0 when array_equal(lb, ub)
# lf1, lf2, uf1, uf2 = r_l2.lb, r_u2.lb, r_l2.ub, r_u2.ub
# Lf1, Lf2, Uf1, Uf2 = Copy(r0.lb), Copy(r0.lb), Copy(r0.ub), Copy(r0.ub)
# Lf1[:, j], Lf2[:, j], Uf1[:, j], Uf2[:, j] = lf1, lf2, uf1, uf2
# r_l2.lb, r_u2.lb, r_l2.ub, r_u2.ub = Lf1, Lf2, Uf1, Uf2
# # if type(r0.definiteRange) not in (bool, bool_) and r0.definiteRange.size != 1:
# # d1, d2 = r_l.definiteRange, r_u.definiteRange
# # D1, D2 = Copy(r0.definiteRange), Copy(r0.definiteRange)
# # D1[:, j], D2[:, j] = d1, d2
# # r_l.definiteRange, r_u.definiteRange = D1, D2
# from numpy.linalg import norm
# if norm(r_u2.ub- r_u.ub) > 1e-7:
# print('uu')
# print(r_u2.ub, r_u.ub)
# print(r_u2.ub - r_u.ub)
## raise 0
# input()
# domain = domain2
## #debug end
for i, v in enumerate(Dep):
domain.modificationVar = v
r_l, r_u = self._iqg(domain, dtype, r0)
if useSlicing and r_l is not r0:# r_l is r0 when array_equal(lb, ub)
lf1, lf2, uf1, uf2 = r_l.lb, r_u.lb, r_l.ub, r_u.ub
Lf1, Lf2, Uf1, Uf2 = Copy(r0.lb), Copy(r0.lb), Copy(r0.ub), Copy(r0.ub)
Lf1[:, j], Lf2[:, j], Uf1[:, j], Uf2[:, j] = lf1, lf2, uf1, uf2
r_l.lb, r_u.lb, r_l.ub, r_u.ub = Lf1, Lf2, Uf1, Uf2
if type(r0.definiteRange) not in (bool, bool_):
print('!')
d1, d2 = r_l.definiteRange, r_u.definiteRange
D1, D2 = atleast_1d(r0.definiteRange).copy(), atleast_1d(r0.definiteRange).copy()
D1[j], D2[j] = d1, d2
r_l.definiteRange, r_u.definiteRange = D1, D2
r[v] = r_l, r_u
if not self.isUncycled:
lf1, lf2, uf1, uf2 = r_l.lb, r_u.lb, r_l.ub, r_u.ub
lf, uf = nanmin(vstack((lf1, lf2)), 0), nanmax(vstack((uf1, uf2)), 0)
if i == 0:
L, U = lf.copy(), uf.copy()
else:
L[Luf] = uf[U>uf].copy()
if not self.isUncycled:
for R in r.values():
r1, r2 = R
if type(r1.lb) != ndarray:
r1.lb, r2.lb, r1.ub, r2.ub = atleast_1d(r1.lb), atleast_1d(r2.lb), atleast_1d(r1.ub), atleast_1d(r2.ub)
r1.lb[r1.lb < L] = L[r1.lb < L]
r2.lb[r2.lb < L] = L[r2.lb < L]
r1.ub[r1.ub > U] = U[r1.ub > U]
r2.ub[r2.ub > U] = U[r2.ub > U]
r0.lb[r0.lb < L] = L[r0.lb < L]
r0.ub[r0.ub > U] = U[r0.ub > U]
# for more safety
domain.useSave = True
domain.useAsMutable = False
domain.modificationVar = None
domain.storedIntervals = {}
return r, r0
def _iqg(self, domain, dtype, r0):
#dep = self._getDep()
v = domain.modificationVar
v_0 = domain[v]
lb, ub = v_0[0], v_0[1]
if v.domain is not None and array_equal(lb, ub):
return r0,r0
assert dtype in (float, float64, int32, int16), 'other types unimplemented yet'
middle = 0.5 * (lb+ub)
if v.domain is not None:
middle1, middle2 = middle.copy(), middle.copy()
adjust_ux_WithDiscreteDomain(middle1, v)
adjust_lx_WithDiscreteDomain(middle2, v)
else:
middle1 = middle2 = middle
domain[v] = (v_0[0], middle1)
domain.localStoredIntervals = {}
r_l = self.interval(domain, dtype, resetStoredIntervals = False)
domain[v] = (middle2, v_0[1])
domain.localStoredIntervals = {}
r_u = self.interval(domain, dtype, resetStoredIntervals = False)
domain[v] = v_0
domain.localStoredIntervals = {}
return r_l, r_u
__pos__ = lambda self: self
# overload "a+b"
# @checkSizes
def __add__(self, other):
# if isinstance(other, Stochastic):
# return other.__add__(self)
for frame_tuple in inspect.stack():
frame = frame_tuple[0]
if 'func_code' in dir(frame) and 'func_code' in dir(npSum) and frame.f_code is npSum.func_code:
pWarn('''
seems like you use numpy.sum() on FuncDesigner object(s),
using FuncDesigner.sum() instead is highly recommended''')
#
## elif frame.f_code is PythonSum.func_code:
## print("python sum!")
#
##
## for S in stk:
## if not S[0].endswith('ooFun.py') and not S[0].endswith('overloads.py') and \
## S[2] == '' and S[3] is not None and 'sum(' in S[3]:
## pWarn('''
## seems like you use Python sum() on FuncDesigner object(s),
## using FuncDesigner.sum() instead is highly recommended''')
if not isinstance(other, (oofun, list, ndarray, tuple)) and not isscalar(other):
raise FuncDesignerException('operation oofun_add is not implemented for the type ' + str(type(other)))
other_is_sum = isinstance(other, oofun) and other._isSum
from overloads import sum
if self._isSum and other_is_sum:
return sum(self._summation_elements + other._summation_elements)
elif self._isSum:
return sum(self._summation_elements + [other])
elif other_is_sum:
return sum(other._summation_elements + [self])
# TODO: check for correct sizes during f, not only f.d
def aux_d(x, y):
Xsize, Ysize = Len(x), Len(y)
if Xsize == 1:
return ones(Ysize)
elif Ysize == 1:
return Eye(Xsize)
elif Xsize == Ysize:
return Eye(Ysize) if not isinstance(x, multiarray) else ones(Ysize).view(multiarray)
else:
raise FuncDesignerException('for oofun summation a+b should be size(a)=size(b) or size(a)=1 or size(b)=1')
if isinstance(other, oofun):
r = oofun(lambda x, y: x+y, [self, other], d = (lambda x, y: aux_d(x, y), lambda x, y: aux_d(y, x)), _isSum = True)
r._summation_elements = [self, other]
r.discrete = self.discrete and other.discrete
r.getOrder = lambda *args, **kwargs: max((self.getOrder(*args, **kwargs), other.getOrder(*args, **kwargs)))
# TODO: move it outside the func, to prevent recreation each time
def interval(domain, dtype):
domain1, definiteRange1 = self._interval(domain, dtype)
domain2, definiteRange2 = other._interval(domain, dtype)
#return domain1[0] + domain2[0], domain1[1] + domain2[1]
return domain1 + domain2, logical_and(definiteRange1, definiteRange2)
r._interval_ = interval
else:
if isscalar(other) and other == 0: return self # sometimes triggers from other parts of FD engine
#if isinstance(other, (OOArray, Stochastic)): return other + self
if isinstance(other, OOArray): return other + self
if isinstance(other, ndarray): other = other.copy()
# if other.size == 1:
# #r = oofun(lambda *ARGS, **KWARGS: None, input = self.input)
# r = oofun(lambda a: a+other, self)
# #r._D = lambda *args, **kwargs: self._D(*args, **kwargs)
# r.d = lambda x: aux_d(x, other)
# #assert len(r._getDep())>0
# #r._c = self._c + other
# else:
r = oofun(lambda a: a+other, self, _isSum=True)
r._summation_elements = [self, other]
r.d = lambda x: aux_d(x, other)
r._getFuncCalcEngine = lambda *args, **kwargs: self._getFuncCalcEngine(*args, **kwargs) + other
r.discrete = self.discrete
r.getOrder = self.getOrder
# TODO: move it outside
Other2 = tile(other, (2, 1))
def interval(domain, dtype):
r, definiteRange = self._interval(domain, dtype)
return r + Other2, definiteRange
r._interval_ = interval
if isscalar(other) or asarray(other).size == 1 or ('size' in self.__dict__ and self.size is asarray(other).size):
r._D = lambda *args, **kwargs: self._D(*args, **kwargs)
r.vectorized = True
return r
__radd__ = lambda self, other: self.__add__(other)
# overload "-a"
def __neg__(self):
if self._neg_elem is not None:
return self._neg_elem
if self._isSum:
from overloads import sum as FDsum
return FDsum([-elem for elem in self._summation_elements])
r = oofun(lambda a: -a, self, d = lambda a: -Eye(Len(a)))
r._neg_elem = self
r._getFuncCalcEngine = lambda *args, **kwargs: -self._getFuncCalcEngine(*args, **kwargs)
r.getOrder = self.getOrder
r._D = lambda *args, **kwargs: dict([(key, -value) for key, value in self._D(*args, **kwargs).items()])
r.d = raise_except
r.criticalPoints = False
r.vectorized = True
def _interval(domain, dtype):
r, definiteRange = self._interval(domain, dtype)
assert r.shape[0] == 2
return -flipud(r), definiteRange
#return (-r[1], -r[0])
r._interval_ = _interval
return r
# overload "a-b"
__sub__ = lambda self, other: self + (-asfarray(other).copy()) if type(other) in (list, tuple, ndarray) else self + (-other)
__rsub__ = lambda self, other: other + (-self)
# overload "a/b"
def __div__(self, other):
if isinstance(other, OOArray):
return other.__rdiv__(self)
if isinstance(other, list): other = asarray(other)
if isscalar(other) or type(other) == ndarray:
return self * (1.0 / other) # to make available using _prod_elements
if isinstance(other, oofun):
r = oofun(lambda x, y: x/y, [self, other])
def aux_dx(x, y):
# TODO: handle float128
y = asfarray(y)
Xsize, Ysize = x.size, y.size
if Xsize != 1:
assert Xsize == Ysize or Ysize == 1, 'incorrect size for oofun devision'
if Xsize != 1:
if Ysize == 1:
r = Diag(None, size=Xsize, scalarMultiplier = 1.0/y)
else:
r = Diag(1.0/y)
else:
r = 1.0 / y
return r
def aux_dy(x, y):
# TODO: handle float128
x = asfarray(x)
Xsize, Ysize = Len(x), Len(y)
r = -x / y**2
if Ysize != 1:
assert Xsize == Ysize or Xsize == 1, 'incorrect size for oofun devision'
r = Diag(r)
return r
r.d = (aux_dx, aux_dy)
def getOrder(*args, **kwargs):
order1, order2 = self.getOrder(*args, **kwargs), other.getOrder(*args, **kwargs)
return order1 if order2 == 0 else inf
r.getOrder = getOrder
def interval(domain, dtype):
lb1_ub1, definiteRange1 = self._interval(domain, dtype)
lb1, ub1 = lb1_ub1[0], lb1_ub1[1]
lb2_ub2, definiteRange2 = other._interval(domain, dtype)
lb2, ub2 = lb2_ub2[0], lb2_ub2[1]
lb2, ub2 = asarray(lb2, dtype), asarray(ub2, dtype)
tmp = vstack((lb1/lb2, lb1/ub2, ub1/lb2, ub1/ub2))
r1, r2 = nanmin(tmp, 0), nanmax(tmp, 0)
ind = logical_or(lb1==0.0, ub1==0.0)
r1[atleast_1d(logical_and(ind, r1>0.0))] = 0.0
r2[atleast_1d(logical_and(ind, r2<0.0))] = 0.0
# adjust inf
# ind1_zero_minus = logical_and(lb1<0, ub1>=0)
# ind1_zero_plus = logical_and(lb1<=0, ub1>0)
ind2_zero_minus = logical_and(lb2<0, ub2>=0)
ind2_zero_plus = logical_and(lb2<=0, ub2>0)
r1[atleast_1d(logical_or(logical_and(ind2_zero_minus, ub1>0), logical_and(ind2_zero_plus, lb1<0)))] = -inf
r2[atleast_1d(logical_or(logical_and(ind2_zero_minus, lb1<0), logical_and(ind2_zero_plus, ub1>0)))] = inf
#assert not any(isnan(r1)) and not any(isnan(r2))
#assert all(r1 <= r2)
return vstack((r1, r2)), logical_and(definiteRange1, definiteRange2)
r._interval_ = interval
else:
# TODO: mb remove it?
other = array(other,'float')# TODO: handle float128
r = oofun(lambda a: a/other, self, discrete = self.discrete)# TODO: involve sparsity if possible!
r.getOrder = self.getOrder
r._getFuncCalcEngine = lambda *args, **kwargs: self._getFuncCalcEngine(*args, **kwargs) / other
#r.d = lambda x: 1.0/other if (isscalar(x) or x.size == 1) else Diag(ones(x.size)/other) if other.size > 1 \
#else Diag(None, size=x.size, scalarMultiplier=1.0/other)
r.d = lambda x: 1.0/other if (isscalar(x) or x.size == 1) else Diag(ones(x.size)/other) #if other.size > 1 \
#else Diag(None, size=x.size, scalarMultiplier=1.0/other)
# commented code is unreacheble, see r._D definition below for other.size == 1
r.criticalPoints = False
# if other.size == 1 or 'size' in self.__dict__ and self.size in (1, other.size):
if other.size == 1:
r._D = lambda *args, **kwargs: dict([(key, value/other) for key, value in self._D(*args, **kwargs).items()])
r.d = raise_except
# r.discrete = self.discrete and (?)
#r.isCostly = True
r.vectorized = True
return r
def __rdiv__(self, other):
# without the code it somehow doesn't fork in either Python3 or latest numpy
#if isinstance(other, Stochastic) or (isinstance(other, OOArray) and any([isinstance(elem, oofun) for elem in atleast_1d(other)])):
if isinstance(other, OOArray) and any([isinstance(elem, oofun) for elem in atleast_1d(other)]):
return other.__div__(self)
other = array(other, 'float') # TODO: sparse matrices handling!
r = oofun(lambda x: other/x, self, discrete = self.discrete)
r.d = lambda x: Diag((- other) / x**2)
def interval(domain, dtype):
arg_lb_ub, definiteRange = self._interval(domain, dtype)
arg_infinum, arg_supremum = arg_lb_ub[0], arg_lb_ub[1]
if other.size != 1:
raise FuncDesignerException('this case for interval calculations is unimplemented yet')
r = vstack((other / arg_supremum, other / arg_infinum))
r1, r2 = amin(r, 0), amax(r, 0)
ind_zero_minus = logical_and(arg_infinum<0, arg_supremum>=0)
r1[atleast_1d(logical_and(ind_zero_minus, other>0))] = -inf
r2[atleast_1d(logical_and(ind_zero_minus, other<0))] = inf
ind_zero_plus = logical_and(arg_infinum<=0, arg_supremum>0)
r1[atleast_1d(logical_and(ind_zero_plus, other<0))] = -inf
r2[atleast_1d(logical_and(ind_zero_plus, other>0))] = inf
return vstack((r1, r2)), definiteRange
r._interval_ = interval
#r.isCostly = True
def getOrder(*args, **kwargs):
order = self.getOrder(*args, **kwargs)
return 0 if order == 0 else inf
r.getOrder = getOrder
r.vectorized = True
return r
# overload "a*b"
def __mul__(self, other):
if isinstance(other, OOArray):#if isinstance(other, (OOArray, Stochastic)):
return other.__mul__(self)
isOtherOOFun = isinstance(other, oofun)
if isinstance(other, list): other = asarray(other)
if self._isProd:
if not isOtherOOFun and not isinstance(self._prod_elements[-1], (oofun, OOArray)):
assert len(self._prod_elements) == 2, 'bug in FD kernel'
return self._prod_elements[0] * (other * self._prod_elements[-1])
if isOtherOOFun:
r = oofun(lambda x, y: x*y, [self, other])
r.d = (lambda x, y: mul_aux_d(x, y), lambda x, y: mul_aux_d(y, x))
r.getOrder = lambda *args, **kwargs: self.getOrder(*args, **kwargs) + other.getOrder(*args, **kwargs)
else:
other = other.copy() if isinstance(other, ndarray) else asarray(other)
r = oofun(lambda x: x*other, self, discrete = self.discrete)
r.getOrder = self.getOrder
r._getFuncCalcEngine = lambda *args, **kwargs: other * self._getFuncCalcEngine(*args, **kwargs)
r.criticalPoints = False
if isscalar(other) or asarray(other).size == 1: # other may be array-like
r._D = lambda *args, **kwargs: dict([(key, value * other) for key, value in self._D(*args, **kwargs).items()])
r.d = raise_except
else:
r.d = lambda x: mul_aux_d(x, other)
r._interval_ = lambda *args, **kw: mul_interval(self, other, isOtherOOFun, *args, **kw)
r.vectorized = True
#r.isCostly = True
r._isProd = True
r._prod_elements = [self, other]
return r
__rmul__ = lambda self, other: self.__mul__(other)
def __pow__(self, other):
if isinstance(other, OOArray):#if isinstance(other, (OOArray, Stochastic)):
return other.__rpow__(self)
# if 'size' not in self.__dict__ or (isscalar(self.size) and self.size == 1):
# return ooarray([self ** other[i] for i in range(other.size)])
# elif isscalar(self.size):# and is not 1
# return ooarray([self[i] ** other[i] for i in range(other.size)])
d_x = lambda x, y: \
(y * x ** (y - 1) if (isscalar(x) or x.size == 1 or isinstance(x, multiarray)) else Diag(y * x ** (y - 1))) if y is not 2 else Diag(2 * x)
d_y = lambda x, y: x ** y * log(x) if (isscalar(y) or y.size == 1) and not isinstance(x, multiarray) else Diag(x ** y * log(x))
if not isinstance(other, oofun):
if isscalar(other):
if type(other) == int: # TODO: handle numpy integer types
pass
#other = asarray(other, dtype='float')
else:
other = asarray(other, dtype= type(other))# with same type, mb float128
elif not isinstance(other, ndarray):
# TODO: fix it wrt int32, int64 etc
other = asarray(other, dtype='float' if type(other) == int else type(other)).copy()
f = lambda x: asanyarray(x) ** other
d = lambda x: d_x(x, other)
input = self
def interval(domain, dtype):
lb_ub, definiteRange = self._interval(domain, dtype)
Tmp = lb_ub ** other
#Tmp = vstack((lb**other, ub**other))
t_min, t_max = nanmin(Tmp, 0), nanmax(Tmp, 0)
#lb2, ub2 = other._interval(domain, dtype) if isinstance(other, oofun) else other, other # TODO: improve it
lb, ub = lb_ub[0], lb_ub[1]
ind = lb < 0.0
if any(ind):
isNonInteger = other != asarray(other, int) # TODO: rational numbers?
# TODO: rework it properly, with matrix operations
if any(isNonInteger):
definiteRange = False
ind_nan = logical_and(logical_and(ind, isNonInteger), ub < 0)
if any(ind_nan):
t_max[atleast_1d(ind_nan)] = nan
#1
t_min[atleast_1d(logical_and(ind, logical_and(t_min>0, ub >= 0)))] = 0.0
# #2
# if asarray(other).size == 1:
# IND = not isNonInteger
# else:
# ind2 = logical_not(isNonInteger)
# IND = other[ind2] % 2 == 0
#
# if any(IND):
# t_min[logical_and(IND, atleast_1d(logical_and(lb<0, ub >= 0)))] = 0.0
return vstack((t_min, t_max)), definiteRange
else:
f = lambda x, y: asanyarray(x) ** y
d = (d_x, d_y)
input = [self, other]
def interval(domain, dtype):
# TODO: handle discrete cases
lb1_ub1, definiteRange1 = self._interval(domain, dtype)
lb1, ub1 = lb1_ub1[0], lb1_ub1[1]
lb2_ub2, definiteRange2 = other._interval(domain, dtype)
lb2, ub2 = lb2_ub2[0], lb2_ub2[1]
T = vstack((lb1 ** lb2, lb1** ub2, ub1**lb1, ub1**ub2))
t_min, t_max = nanmin(T, 0), nanmax(T, 0)
definiteRange = logical_and(definiteRange1, definiteRange2)
ind1 = lb1 < 0
# TODO: check it, especially with integer "other"
if any(ind1):
# TODO: rework it with matrix operations
definiteRange = False
ind2 = ub1 >= 0
t_min[atleast_1d(logical_and(logical_and(ind1, ind2), logical_and(t_min > 0.0, ub2 > 0.0)))] = 0.0
t_max[atleast_1d(logical_and(ind1, logical_not(ind2)))] = nan
t_min[atleast_1d(logical_and(ind1, logical_not(ind2)))] = nan
return vstack((t_min, t_max)), definiteRange
r = oofun(f, input, d = d, _interval_=interval)
if isinstance(other, oofun) or (not isinstance(other, int) or (type(other) == ndarray and other.flatten()[0] != int)):
r.attach((self>0)('pow_domain_%d'%r._id, tol=-1e-7)) # TODO: if "other" is fixed oofun with integer value - omit this
r.isCostly = True
r.vectorized = True
return r
def __rpow__(self, other):
assert not isinstance(other, oofun)# if failed - check __pow__implementation
if isscalar(other):
if type(other) == int: # TODO: handle numpy integer types
other = float(other)
elif not isinstance(other, ndarray):
other = asarray(other, 'float' if type(other) in (int, int32, int64, int16, int8) else type(other))
f = lambda x: other ** x
#d = lambda x: Diag(asarray(other) **x * log(asarray(other)))
d = lambda x: Diag(other **x * log(other))
r = oofun(f, self, d=d, criticalPoints = False)
#r.isCostly = True
r.vectorized = True
return r
def __xor__(self, other): raise FuncDesignerException('For power of oofuns use a**b, not a^b')
def __rxor__(self, other): raise FuncDesignerException('For power of oofuns use a**b, not a^b')
def __getitem__(self, ind): # overload for oofun[ind]
if isinstance(ind, oofun):# NOT IMPLEMENTED PROPERLY YET
self.pWarn('Slicing oofun by oofun IS NOT IMPLEMENTED PROPERLY YET')
f = lambda x, _ind: x[_ind]
def d(x, _ind):
r = zeros(x.shape)
r[_ind] = 1
return r
elif type(ind) not in (int, int32, int64, int16, int8):
# Python 3 slice
return self.__getslice__(ind.start, ind.stop)
else:
if not hasattr(self, '_slicesIndexDict'):
self._slicesIndexDict = {}
if ind in self._slicesIndexDict:
return self._slicesIndexDict[ind]
f = lambda x: x[ind]
def d(x):
Xsize = Len(x)
condBigMatrix = Xsize > 100
if condBigMatrix and scipyInstalled:
r = SparseMatrixConstructor((1, x.shape[0]))
r[0, ind] = 1.0
else:
if condBigMatrix and not scipyInstalled: self.pWarn(scipyAbsentMsg)
r = zeros_like(x)
r[ind] = 1
return r
r = oofun(f, self, d = d, size = 1, getOrder = self.getOrder)
# TODO: check me!
# what about a[a.size/2:]?
# TODO: edit me!
# if self.is_oovar:
# r.is_oovarSlice = True
self._slicesIndexDict[ind] = r
return r
def __getslice__(self, ind1, ind2):# overload for oofun[ind1:ind2]
#TODO: mb check if size is known then use it instead of None?
if ind1 is None:
ind1 = 0
if ind2 is None:
if 'size' in self.__dict__ and type(self.size) in (int, int8, int16, int32, int64):
ind2 = self.size
else:
raise FuncDesignerException('if oofun.size is not provided then you should provide full slice coords, e.g. x[3:10], not x[3:]')
assert not isinstance(ind1, oofun) and not isinstance(ind2, oofun), 'slicing by oofuns is unimplemented yet'
f = lambda x: x[ind1:ind2]
def d(x):
condBigMatrix = Len(x) > 100 #and (ind2-ind1) > 0.25*x.size
if condBigMatrix and not scipyInstalled:
self.pWarn(scipyAbsentMsg)
if condBigMatrix and scipyInstalled:
r = SP_eye(ind2-ind1, ind2-ind1)
if ind1 != 0:
m1 = SparseMatrixConstructor((ind2-ind1, ind1))
r = Hstack((SparseMatrixConstructor((ind2-ind1, ind1)), r))
if ind2 != x.size:
r = Hstack((r, SparseMatrixConstructor((ind2-ind1, x.size - ind2))))
else:
m1 = zeros((ind2-ind1, ind1))
m2 = eye(ind2-ind1)
m3 = zeros((ind2-ind1, x.size - ind2))
r = hstack((m1, m2, m3))
return r
r = oofun(f, self, d = d, getOrder = self.getOrder)
return r
#def __len__(self):
#return self.size
#raise FuncDesignerException('using len(obj) (where obj is oovar or oofun) is not possible (at least yet), use obj.size instead')
def sum(self):
def d(x):
if type(x) == ndarray and x.ndim > 1: raise FuncDesignerException('sum(x) is not implemented yet for arrays with ndim > 1')
return ones_like(x)
def interval(domain, dtype):
if type(domain) == ooPoint and domain.isMultiPoint:
raise FuncDesignerException('interval calculations are unimplemented for sum(oofun) yet')
lb_ub, definiteRange = self._interval(domain, dtype)
lb, ub = lb_ub[0], lb_ub[1]
return vstack((npSum(lb, 0), npSum(ub, 0))), definiteRange
r = oofun(npSum, self, getOrder = self.getOrder, _interval_ = interval, d=d)
return r
def prod(self):
# TODO: consider using r.isCostly = True
r = oofun(prod, self)
#r.getOrder = lambda *args, **kwargs: self.getOrder(*args, **kwargs)*self.size
def d(x):
x = asarray(x) # prod is used rarely, so optimizing it is not important
if x.ndim > 1: raise FuncDesignerException('prod(x) is not implemented yet for arrays with ndim > 1')
ind_zero = where(x==0)[0].tolist()
ind_nonzero = nonzero(x)[0].tolist()
numOfZeros = len(ind_zero)
r = prod(x) / x
if numOfZeros >= 2:
r[ind_zero] = 0
elif numOfZeros == 1:
r[ind_zero] = prod(x[ind_nonzero])
return r
r.d = d
return r
""" Handling constraints """
# TODO: optimize for lb-ub imposed on oovars
# TODO: fix it for discrete problems like MILP, MINLP
def __gt__(self, other): # overload for >
if self.is_oovar and not isinstance(other, (oofun, OOArray)) and not (isinstance(other, ndarray) and str(other.dtype) =='object'):
r = BoxBoundConstraint(self, lb = other)
elif isinstance(other, OOArray) or (isinstance(other, ndarray) and str(other.dtype) =='object'):
return other.__le__(self)
else:
r = Constraint(self - other, lb=0.0) # do not perform check for other == 0, copy should be returned, not self!
return r
# overload for >=
__ge__ = __gt__
# TODO: fix it for discrete problems like MILP
def __lt__(self, other): # overload for <
# TODO:
#(self.is_oovar or self.is_oovarSlice)
if self.is_oovar and not isinstance(other, (oofun, OOArray)) and not(isinstance(other, ndarray) and str(other.dtype) =='object'):
r = BoxBoundConstraint(self, ub = other)
elif isinstance(other, OOArray) or (isinstance(other, ndarray) and str(other.dtype) =='object'):
return other.__ge__(self)
else:
r = Constraint(self - other, ub = 0.0) # do not perform check for other == 0, copy should be returned, not self!
return r
# overload for <=
__le__ = __lt__
__eq__ = lambda self, other: self.eq(other)
def eq(self, other):
#if other in (None, (), []): return False
if other is None or other is () or (type(other) == list and len(other) == 0): return False
if type(other) in (str, string_):
#if self.domain is not None and self.domain is not bool and self.domain is not 'bool':
if 'aux_domain' not in self.__dict__:
if not self.is_oovar:
raise FuncDesignerException('comparing with non-numeric data is allowed for string oovars, not for oofuns')
self.formAuxDomain()
# if len(self.domain) != len(self.aux_domain):
# raise FuncDesignerException('probably you have changed domain of categorical oovar, that is not allowed')
ind = searchsorted(self.aux_domain, other, 'left')
if self.aux_domain[ind] != other:
raise FuncDesignerException('compared value %s is absent in oovar %s domain' %(other, self.name))
#r = (self == ind)(tol=0.5)
r = Constraint(self - ind, ub = 0.0, lb = 0.0, tol=0.5)
#print (self, other)
# TODO: use it instead
if self.is_oovar: r.nlh = lambda Lx, Ux, p, dataType: self.nlh(Lx, Ux, p, dataType, ind)
return r
if 'startswith' in dir(other): return False
#if self.is_oovar and not isinstance(other, oofun):
#raise FuncDesignerException('Constraints like this: "myOOVar = " are not implemented yet and are not recommended; for openopt use freeVars / fixedVars instead')
r = Constraint(self - other, ub = 0.0, lb = 0.0) # do not perform check for other == 0, copy should be returned, not self!
if self.is_oovar and isscalar(other) and self.domain is not None:
if self.domain is bool or self.domain is 'bool':
if other not in [0, 1]:# and type(other) not in (int, int16, int32, int64):
raise FuncDesignerException('bool oovar can be compared with [0,1] only')
r.nlh = self.nlh if other == 1.0 else (~self).nlh
r.alt_nlh_func = True
elif self.domain is not int and self.domain is not 'int':# and type(other) in (str, string_):
pass
# if self.is_oovar:
# self.formAuxDomain()
# print('2', other)
# ind = searchsorted(self.aux_domain, other, 'left')
# #r.nlh = lambda Lx, Ux, p, dataType: self.nlh(Lx, Ux, p, dataType, ind)
# r = Constraint(self - ind, ub = 0.0, lb = 0.0)
# r.nlh = lambda Lx, Ux, p, dataType: self.nlh(Lx, Ux, p, dataType, ind)
# print (self.domain, other)
# r_nlh = r.nlh
# def nlh_c(Lx, Ux, p, dataType):
# r1 = nlh2(Lx, Ux, p, dataType)
# r2 = r_nlh(Lx, Ux, p, dataType)
# print(r1)
# print(r2)
# raw_input()
# #r.nlh = lambda Lx, Ux, p, dataType: self.nlh(Lx, Ux, p, dataType, other)
# r.nlh = nlh_c
## r.nlh = lambda Lx, Ux, p, dataType: self.nlh(Lx, Ux, p, dataType, other)
## r.alt_nlh_func = True
return r
""" getInput """
def _getInput(self, *args, **kwargs):
# self.inputOOVarTotalLength = 0
r = []
for item in self.input:
tmp = item._getFuncCalcEngine(*args, **kwargs) if isinstance(item, oofun) else item(*args, **kwargs) if isinstance(item, OOArray) else item
r.append(tmp if type(tmp) not in (list, tuple, Stochastic) else asanyarray(tmp))
return tuple(r)
""" getDep """
def _getDep(self):
# returns Python set of oovars it depends on
if hasattr(self, 'dep'):
return self.dep
elif self.input is None:
self.dep = None
else:
if type(self.input) not in (list, tuple) and not isinstance(self.input, OOArray):
self.input = [self.input]
#OLD
# r = set()
# for oofunInstance in self.input:
# if not isinstance(oofunInstance, oofun): continue
# if oofunInstance.is_oovar:
# r.add(oofunInstance)
# continue
# tmp = oofunInstance._getDep()
# if tmp is None: continue
# r.update(tmp)
# self.dep = r
# / OLD
# NEW
r_oovars = []
r_oofuns = []
isUncycled = True
Tmp = set()
for Elem in self.input:
if isinstance(Elem, OOArray):
for _elem in Elem:
if isinstance(_elem, oofun):
Tmp.add(_elem)
# _tmp = _elem._getDep()
# if _tmp is not None or (isinstance(_tmp, oofun) and _tmp.is_oovar):
# Tmp.add(_tmp)
#Tmp.update([[_elem._getDep() for _elem in Elem] for Elem in self.input if isinstance(Elem, OOArray)])
for Elem in (list(Tmp) + self.input):
if not isinstance(Elem, oofun): continue
if Elem.is_oovar:
r_oovars.append(Elem)
continue
tmp = Elem._getDep()
if not Elem.isUncycled: isUncycled = False
if tmp is None or len(tmp)==0: continue # TODO: remove None, use [] instead
r_oofuns.append(tmp)
r = set(r_oovars)
# Python 2.5 sel.update fails on empty input
if len(r_oofuns)!=0: r.update(*r_oofuns)
if len(r_oovars) + sum([len(elem) for elem in r_oofuns]) != len(r):
isUncycled = False
self.isUncycled = isUncycled
self.dep = r
# /NEW
return self.dep
""" getFunc """
def _getFunc(self, *args, **kwargs):
Args = args
if len(args) == 0 and len(kwargs) == 0:
raise FuncDesignerException('at least one argument is required')
if len(args) != 0:
if type(args[0]) != str:
assert not isinstance(args[0], oofun), "you can't invoke oofun on another one oofun"
x = args[0]
if isinstance(x, dict) and not isinstance(x, ooPoint):
x = ooPoint(x)
Args = (x,)+args[1:]
if self.is_oovar:
return self._getFuncCalcEngine(*Args, **kwargs)
# if isinstance(x, dict):
# tmp = x.get(self, None)
# if tmp is not None:
# # currently tmp hasn't to be sparse matrix, mb for future
# if isinstance(tmp, Stochastic):
# r = getattr(x, 'maxDistributionSize', inf)
# tmp.maxDistributionSize = r
# return tmp
# else:
# return float(tmp) if isscalar(tmp) and type(tmp)==int else asfarray(tmp) if not isspmatrix(tmp) else tmp
# elif self.name in x:
# tmp = x[self.name]
# return float(tmp) if isscalar(tmp) and type(tmp)==int else asfarray(tmp) if not isspmatrix(tmp) else tmp
# else:
# s = 'for oovar ' + self.name + \
# " the point involved doesn't contain neither name nor the oovar instance. Maybe you try to get function value or derivative in a point where value for an oovar is missing"
# raise FuncDesignerException(s)
# elif hasattr(x, 'xf'):
# if x.probType == 'MOP': # x is MOP result struct
# s = 'evaluation of MOP result on arguments is unimplemented yet, use r.solutions'
# raise FuncDesignerException(s)
# # TODO: possibility of squeezing
# return x.xf[self]
# else:
# raise FuncDesignerException('Incorrect data type (%s) while obtaining oovar %s value' %(type(x), self.name))
else:
self.name = args[0]
return self
else:
for fn in ['name', 'size', 'tol']:
if fn in kwargs:
setattr(self, fn, kwargs[fn])
return self
if hasattr(x, 'probType') and x.probType == 'MOP':# x is MOP result struct
s = 'evaluation of MOP result on arguments is unimplemented yet, use r.solutions'
raise FuncDesignerException(s)
return self._getFuncCalcEngine(*Args, **kwargs)
def _getFuncCalcEngine(self, *args, **kwargs):
x = args[0]
dep = self._getDep()
CondSamePointByID = True if type(x) == ooPoint and not x.isMultiPoint and self._point_id == x._id else False
fixedVarsScheduleID = kwargs.get('fixedVarsScheduleID', -1)
fixedVars = kwargs.get('fixedVars', None)
Vars = kwargs.get('Vars', None)
sameVarsScheduleID = fixedVarsScheduleID == self._lastFuncVarsID
rebuildFixedCheck = not sameVarsScheduleID
if fixedVarsScheduleID != -1: self._lastFuncVarsID = fixedVarsScheduleID
if rebuildFixedCheck:
self._isFixed = (fixedVars is not None and dep.issubset(fixedVars)) or (Vars is not None and dep.isdisjoint(Vars))
if isinstance(x, ooPoint) and x.isMultiPoint:
cond_same_point = False
else:
cond_same_point = CondSamePointByID or \
(self._f_val_prev is not None and (self._isFixed or (self.isCostly and all([array_equal((x if isinstance(x, dict) else x.xf)[elem], self._f_key_prev[elem]) for elem in (dep & set((x if isinstance(x, dict) else x.xf).keys()))]))))
if cond_same_point:
self.same += 1
tmp = self._f_val_prev
return tmp.copy() if isinstance(tmp, (ndarray, Stochastic)) else tmp
self.evals += 1
#TODO: add condition "and self in x._p.dictOfLinearFuncs" instead of self._order == 1
#use_line_points = hasattr(x,'_p') and x._p.solver.useLinePoints and self._order == 1
# if use_line_points:
# _linePointDescriptor = getattr(x, '_linePointDescriptor', None)
# if _linePointDescriptor is not None:
# #point1, alp, point2 = _linePointDescriptor
# alp = _linePointDescriptor
# r1, r2 = self._p._firstLinePointDict[self], self._p._secondLinePointDict[self]
# #assert r1 is not None and r2 is not None
# return r1 * (1-alp) + r2 * alp
if type(self.args) != tuple:
self.args = (self.args, )
Input = self._getInput(*args, **kwargs)
# if not isinstance(x, ooPoint) or not x.isMultiPoint or (self.vectorized and not any([isinstance(inp, Stochastic) for inp in Input])):
if not isinstance(x, ooPoint) or not x.isMultiPoint or self.vectorized:
if self.args != ():
Input += self.args
Tmp = self.fun(*Input)
if isinstance(Tmp, (list, tuple)):
tmp = hstack(Tmp) if len(Tmp) > 1 else Tmp[0]
else:
tmp = Tmp
else:
if hasattr(x, 'N'):
N = x.N
else:
# TODO: fix it for x.values() is Stochastic
N = 1
for inp in Input:
if not isinstance(inp, Stochastic):
N = inp.size if type(inp) == ndarray else 1
break
inputs = zip(*[(atleast_1d(inp) if not isinstance(inp, Stochastic) and type(inp) == ndarray and inp.size == N else [inp]*N) for inp in Input])
# Check it!
Tmp = [self.fun(*inp) if self.args == () else self.fun(*(inp + self.args)) for inp in inputs]
if N == 1:
tmp = Tmp[0]
else:
tmp = array([elem for elem in Tmp], object).view(multiarray)
#if self._c != 0.0: tmp += self._c
#self.outputTotalLength = ([asarray(elem).size for elem in self.fun(*Input)])#self.f_val_prev.size # TODO: omit reassigning
#!! TODO: handle case tmp is multiarray of Stochastic
if isinstance(tmp, Stochastic):
if 'xf' in x.__dict__:
maxDistributionSize = getattr(x.xf, 'maxDistributionSize', 0)
else:
maxDistributionSize = getattr(x, 'maxDistributionSize', 0)
if maxDistributionSize == 0:
s = '''
if one of function arguments is stochastic distribution
without resolving into quantified value
(e.g. uniform(-10,10) instead of uniform(-10,10, 100), 100 is number of point to emulate)
then you should evaluate the function
onto oopoint with assigned parameter maxDistributionSize'''
raise FuncDesignerException(s)
if tmp.size > maxDistributionSize:
tmp.reduce(maxDistributionSize)
tmp.maxDistributionSize = maxDistributionSize
if ((type(x) == ooPoint and not x.isMultiPoint) and not (isinstance(tmp, ndarray) and type(tmp) != ndarray)) or self._isFixed:# or self.isCostly:
# TODO: rework it (for input with ooarays)
try:
self._f_key_prev = dict([(elem, copy((x if isinstance(x, dict) else x.xf)[elem])) for elem in dep]) if self.isCostly else None
self._f_val_prev = tmp.copy() if isinstance(tmp, (ndarray, Stochastic)) else tmp
if type(x) == ooPoint:
self._point_id = x._id
except:
pass
r = tmp
# if use_line_points:
# self._p._currLinePointDict[self] = r
# if fixedVarsScheduleID != -1:
# self._lastSize = tmp.size
return r
""" getFunc """
__call__ = lambda self, *args, **kwargs: self._getFunc(*args, **kwargs)
""" derivatives """
def D(self, x, Vars=None, fixedVars = None, resultKeysType = 'vars', useSparse = False, exactShape = False, fixedVarsScheduleID = -1):
# resultKeysType doesn't matter for the case isinstance(Vars, oovar)
if Vars is not None and fixedVars is not None:
raise FuncDesignerException('No more than one argument from "Vars" and "fixedVars" is allowed for the function')
#assert type(Vars) != ndarray and type(fixedVars) != ndarray
if not isinstance(x, ooPoint): x = ooPoint(x)
initialVars = Vars
#TODO: remove cloned code
if Vars is not None:
if type(Vars) in [list, tuple]:
Vars = set(Vars)
elif isinstance(Vars, oofun):
if not Vars.is_oovar:
raise FuncDesignerException('argument Vars is expected as oovar or python list/tuple of oovar instances')
Vars = set([Vars])
if fixedVars is not None:
if type(fixedVars) in [list, tuple]:
fixedVars = set(fixedVars)
elif isinstance(fixedVars, oofun):
if not fixedVars.is_oovar:
raise FuncDesignerException('argument fixedVars is expected as oovar or python list/tuple of oovar instances')
fixedVars = set([fixedVars])
r = self._D(x, fixedVarsScheduleID, Vars, fixedVars, useSparse = useSparse)
r = dict([(key, (val if type(val)!=DiagonalType else val.resolve(useSparse))) for key, val in r.items()])
is_oofun = isinstance(initialVars, oofun)
if is_oofun and not initialVars.is_oovar:
# TODO: handle it with input of type list/tuple/etc as well
raise FuncDesignerException('Cannot perform differentiation by non-oovar input')
if resultKeysType == 'names':
raise FuncDesignerException("""This possibility is out of date,
if it is still present somewhere in FuncDesigner doc inform developers""")
elif resultKeysType == 'vars':
rr = {}
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! TODO: remove the cycle!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
for oov, val in x.items():
#if not isinstance(oov not in r, bool): print oov not in r
if oov not in r or (fixedVars is not None and oov in fixedVars):
continue
tmp = r[oov]
if useSparse == False and hasattr(tmp, 'toarray'): tmp = tmp.toarray()
if not exactShape and not isspmatrix(tmp) and not isscalar(tmp):
if tmp.size == 1: tmp = asscalar(tmp)
elif min(tmp.shape) == 1: tmp = tmp.flatten()
rr[oov] = tmp
return rr if not is_oofun else rr[initialVars]
else:
raise FuncDesignerException('Incorrect argument resultKeysType, should be "vars" or "names"')
def _D(self, x, fixedVarsScheduleID, Vars=None, fixedVars = None, useSparse = 'auto'):
if self.is_oovar:
if (fixedVars is not None and self in fixedVars) or (Vars is not None and self not in Vars):
return {}
tmp = x[self]
return {self:Eye(asarray(tmp).size)} if not isinstance(tmp, multiarray) else {self: ones_like(tmp).view(multiarray)}
if self.input[0] is None: return {} # fixed oofun. TODO: implement input = [] properly
if self.discrete:
return {}
#raise FuncDesignerException('The oofun or oovar instance has been declared as discrete, no derivative is available')
CondSamePointByID = True if isinstance(x, ooPoint) and self._point_id1 == x._id else False
sameVarsScheduleID = fixedVarsScheduleID == self._lastDiffVarsID
dep = self._getDep()
rebuildFixedCheck = not sameVarsScheduleID
if rebuildFixedCheck:
self._isFixed = (fixedVars is not None and dep.issubset(fixedVars)) or (Vars is not None and dep.isdisjoint(Vars))
if self._isFixed: return {}
##########################
# TODO: optimize it. Omit it for simple cases.
#isTransmit = self._usedIn == 1 # Exactly 1! not 0, 2, ,3, 4, etc
#involveStore = not isTransmit or self._directlyDwasInwolved
involveStore = self.isCostly
#cond_same_point = hasattr(self, '_d_key_prev') and sameDerivativeVariables and (CondSamePointByID or (involveStore and all([array_equal(x[elem], self.d_key_prev[elem]) for elem in dep])))
cond_same_point = sameVarsScheduleID and \
((CondSamePointByID and self._d_val_prev is not None) or \
(involveStore and self._d_key_prev is not None and all([array_equal(x[elem], self._d_key_prev[elem]) for elem in dep])))
if cond_same_point:
self.same_d += 1
#return deepcopy(self.d_val_prev)
return dict([(key, Copy(val)) for key, val in self._d_val_prev.items()])
else:
self.evals_d += 1
if isinstance(x, ooPoint): self._point_id1 = x._id
if fixedVarsScheduleID != -1: self._lastDiffVarsID = fixedVarsScheduleID
derivativeSelf = self._getDerivativeSelf(x, fixedVarsScheduleID, Vars, fixedVars)
r = Derivative()
ac = -1
for i, inp in enumerate(self.input):
if not isinstance(inp, oofun): continue
if inp.discrete: continue
if inp.is_oovar:
if (Vars is not None and inp not in Vars) or (fixedVars is not None and inp in fixedVars):
continue
ac += 1
tmp = derivativeSelf[ac]
val = r.get(inp, None)
if val is not None:
if isscalar(tmp) or (type(val) == type(tmp) == ndarray and prod(tmp.shape) <= prod(val.shape)): # some sparse matrices has no += implemented
r[inp] += tmp
else:
if isspmatrix(val) and type(tmp) == DiagonalType:
tmp = tmp.resolve(True)
r[inp] = r[inp] + tmp
else:
r[inp] = tmp
else:
ac += 1
elem_d = inp._D(x, fixedVarsScheduleID, Vars=Vars, fixedVars=fixedVars, useSparse = useSparse)
t1 = derivativeSelf[ac]
for key, val in elem_d.items():
#if isinstance(t1, Stochastic) or isinstance(val, Stochastic):
#rr = t1 * val
if isinstance(t1, Stochastic) or ((isscalar(val) or isinstance(val, multiarray)) and (isscalar(t1) or isinstance(t1, multiarray))):
rr = t1 * val
elif isinstance(val, Stochastic):
rr = val * t1
elif type(t1) == DiagonalType and type(val) == DiagonalType:
rr = t1 * val
elif type(t1) == DiagonalType or type(val) == DiagonalType:
if isspmatrix(t1): # thus val is DiagonalType
rr = t1._mul_sparse_matrix(val.resolve(True))
else:
if not isPyPy or type(val) != DiagonalType:
rr = t1 * val #if type(t1) == DiagonalType or type(val) not in (ndarray, DiagonalType) else (val.T * t1).T # for PyPy compatibility
else:
rr = (val * t1.T).T
elif isscalar(val) or isscalar(t1) or prod(t1.shape)==1 or prod(val.shape)==1:
rr = (t1 if isscalar(t1) or prod(t1.shape)>1 else asscalar(t1) if isinstance(t1, ndarray) else t1[0, 0]) \
* (val if isscalar(val) or prod(val.shape)>1 else asscalar(val) if isinstance(val, ndarray) else val[0, 0])
else:
if val.ndim < 2: val = atleast_2d(val)
if useSparse is False:
t2 = val
else:
t1, t2 = self._considerSparse(t1, val)
if not type(t1) == type(t2) == ndarray:
# CHECKME: is it trigger somewhere?
if not scipyInstalled:
self.pWarn(scipyAbsentMsg)
rr = dot(t1, t2)
else:
t1 = t1 if isspmatrix_csc(t1) else t1.tocsc() if isspmatrix(t1) else scipy.sparse.csc_matrix(t1)
t2 = t2 if isspmatrix_csr(t2) else t2.tocsr() if isspmatrix(t2) else scipy.sparse.csr_matrix(t2)
if t2.shape[0] != t1.shape[1]:
if t2.shape[1] == t1.shape[1]:
t2 = t2.T
else:
raise FuncDesignerException('incorrect shape in FuncDesigner function _D(), inform developers about the bug')
rr = t1._mul_sparse_matrix(t2)
if useSparse is False:
rr = rr.toarray()
else:
rr = dot(t1, t2)
#assert rr.ndim>1
Val = r.get(key, None)
ValType = type(Val)
if Val is not None:
if type(rr) == DiagonalType:
if ValType == DiagonalType:
Val = Val + rr # !!!! NOT inplace! (elseware will overwrite stored data used somewhere else)
else:
tmp = rr.resolve(useSparse)
if type(tmp) == ndarray and hasattr(Val, 'toarray'):
Val = Val.toarray()
if type(tmp) == ValType == ndarray and Val.size >= tmp.size:
Val += tmp
else: # may be problems with sparse matrices inline operation, which are badly done in scipy.sparse for now
Val = Val + tmp
else:
if isinstance(Val, ndarray) and hasattr(rr, 'toarray'): # i.e. rr is sparse matrix
rr = rr.toarray() # I guess r[key] will hardly be all-zeros
elif hasattr(Val, 'toarray') and isinstance(rr, ndarray): # i.e. Val is sparse matrix
Val = Val.toarray()
if type(rr) == ValType == ndarray and rr.size == Val.size:
Val += rr
else:
Val = Val + rr
r[key] = Val
else:
r[key] = rr
self._d_val_prev = dict([(key, Copy(value)) for key, value in r.items()])
self._d_key_prev = dict([(elem, Copy(x[elem])) for elem in dep]) if involveStore else None
return r
# TODO: handle 2**15 & 0.25 as parameters
def _considerSparse(self, t1, t2):
if int64(prod(t1.shape)) * int64(prod(t2.shape)) > 2**15 and (isinstance(t1, ndarray) and t1.nonzero()[0].size < 0.25*t1.size) or \
(isinstance(t2, ndarray) and t2.nonzero()[0].size < 0.25*t2.size):
if scipy is None:
self.pWarn(scipyAbsentMsg)
return t1, t2
if not isinstance(t1, scipy.sparse.csc_matrix):
t1 = scipy.sparse.csc_matrix(t1)
if t1.shape[1] != t2.shape[0]: # can be from flattered t1
assert t1.shape[0] == t2.shape[0], 'bug in FuncDesigner Kernel, inform developers'
t1 = t1.T
if not isinstance(t2, scipy.sparse.csr_matrix):
t2 = scipy.sparse.csr_matrix(t2)
return t1, t2
def _getDerivativeSelf(self, x, fixedVarsScheduleID, Vars, fixedVars):
Input = self._getInput(x, fixedVarsScheduleID=fixedVarsScheduleID, Vars=Vars, fixedVars=fixedVars)
expectedTotalInputLength = sum([Len(elem) for elem in Input])
# if hasattr(self, 'size') and isscalar(self.size): nOutput = self.size
# else: nOutput = self(x).size
hasUserSuppliedDerivative = self.d is not None
if hasUserSuppliedDerivative:
derivativeSelf = []
if type(self.d) == tuple:
if len(self.d) != len(self.input):
raise FuncDesignerException('oofun error: num(derivatives) not equal to neither 1 nor num(inputs)')
for i, deriv in enumerate(self.d):
inp = self.input[i]
if not isinstance(inp, oofun) or inp.discrete:
#if deriv is not None:
#raise FuncDesignerException('For an oofun with some input oofuns declared as discrete you have to set oofun.d[i] = None')
continue
#!!!!!!!!! TODO: handle fixed cases properly!!!!!!!!!!!!
#if hasattr(inp, 'fixed') and inp.fixed: continue
if inp.is_oovar and ((Vars is not None and inp not in Vars) or (fixedVars is not None and inp in fixedVars)):
continue
if deriv is None:
if not DerApproximatorIsInstalled:
raise FuncDesignerException('To perform gradients check you should have DerApproximator installed, see http://openopt.org/DerApproximator')
derivativeSelf.append(get_d1(self.fun, Input, diffInt=self.diffInt, stencil = self.stencil, \
args=self.args, varForDifferentiation = i, pointVal = self._getFuncCalcEngine(x), exactShape = True))
else:
# !!!!!!!!!!!!!! TODO: add check for user-supplied derivative shape
tmp = deriv(*Input)
if not isscalar(tmp) and type(tmp) in (ndarray, tuple, list) and type(tmp) != DiagonalType: # i.e. not a scipy.sparse matrix
tmp = atleast_2d(tmp)
########################################
_tmp = Input[i]
Tmp = 1 if isscalar(_tmp) or prod(_tmp.shape) == 1 else len(Input[i])
if tmp.shape[1] != Tmp:
# TODO: add debug msg
# print('incorrect shape in FD AD _getDerivativeSelf')
# print tmp.shape[0], nOutput, tmp
if tmp.shape[0] != Tmp: raise FuncDesignerException('error in getDerivativeSelf()')
tmp = tmp.T
########################################
derivativeSelf.append(tmp)
else:
tmp = self.d(*Input)
if not isscalar(tmp) and type(tmp) in (ndarray, tuple, list): # i.e. not a scipy.sparse matrix
tmp = atleast_2d(tmp)
if tmp.shape[1] != expectedTotalInputLength:
# TODO: add debug msg
if tmp.shape[0] != expectedTotalInputLength: raise FuncDesignerException('error in getDerivativeSelf()')
tmp = tmp.T
ac = 0
if isinstance(tmp, ndarray) and hasattr(tmp, 'toarray') and not isinstance(tmp, multiarray): tmp = tmp.A # is dense matrix
#if not isinstance(tmp, ndarray) and not isscalar(tmp) and type(tmp) != DiagonalType:
if len(Input) == 1:
# if type(tmp) == DiagonalType:
# # TODO: mb rework it
# if Input[0].size > 150 and tmp.size > 150:
# tmp = tmp.resolve(True).tocsc()
# else: tmp = tmp.resolve(False)
derivativeSelf = [tmp]
else:
for i, inp in enumerate(Input):
t = self.input[i]
if t.discrete or (t.is_oovar and ((Vars is not None and t not in Vars) or (fixedVars is not None and t in fixedVars))):
ac += inp.size
continue
if isinstance(tmp, ndarray):
TMP = tmp[:, ac:ac+Len(inp)]
elif isscalar(tmp):
TMP = tmp
elif type(tmp) == DiagonalType:
if tmp.size == inp.size and ac == 0:
TMP = tmp
else:
# print debug warning here
# TODO: mb rework it
if inp.size > 150 and tmp.size > 150:
tmp = tmp.resolve(True).tocsc()
else: tmp = tmp.resolve(False)
TMP = tmp[:, ac:ac+inp.size]
else: # scipy.sparse matrix
TMP = tmp.tocsc()[:, ac:ac+inp.size]
ac += Len(inp)
derivativeSelf.append(TMP)
# TODO: is it required?
# if not hasattr(self, 'outputTotalLength'): self(x)
#
# if derivativeSelf.shape != (self.outputTotalLength, self.inputTotalLength):
# s = 'incorrect shape for user-supplied derivative of oofun '+self.name+': '
# s += '(%d, %d) expected, (%d, %d) obtained' % (self.outputTotalLength, self.inputTotalLength, derivativeSelf.shape[0], derivativeSelf.shape[1])
# raise FuncDesignerException(s)
else:
if Vars is not None or fixedVars is not None: raise FuncDesignerException("sorry, custom oofun derivatives don't work with Vars/fixedVars arguments yet")
if not DerApproximatorIsInstalled:
raise FuncDesignerException('To perform this operation you should have DerApproximator installed, see http://openopt.org/DerApproximator')
derivativeSelf = get_d1(self.fun, Input, diffInt=self.diffInt, stencil = self.stencil, args=self.args, pointVal = self._getFuncCalcEngine(x), exactShape = True)
if type(derivativeSelf) == tuple:
derivativeSelf = list(derivativeSelf)
elif type(derivativeSelf) != list:
derivativeSelf = [derivativeSelf]
#assert all([elem.ndim > 1 for elem in derivativeSelf])
# assert len(derivativeSelf[0])!=16
#assert (type(derivativeSelf[0]) in (int, float)) or derivativeSelf[0][0]>480.00006752 or derivativeSelf[0][0]<480.00006750
return derivativeSelf
def D2(self, x):
raise FuncDesignerException('2nd derivatives for obj-funcs are not implemented yet')
def check_d1(self, point):
if self.d is None:
self.disp('Error: no user-provided derivative(s) for oofun ' + self.name + ' are attached')
return # TODO: return non-void result
separator = 75 * '*'
self.disp(separator)
assert type(self.d) != list
val = self(point)
input = self._getInput(point)
ds= self._getDerivativeSelf(point, fixedVarsScheduleID = -1, Vars=None, fixedVars=None)
self.disp(self.name + ': checking user-supplied gradient')
self.disp('according to:')
self.disp(' diffInt = ' + str(self.diffInt)) # TODO: ADD other parameters: allowed epsilon, maxDiffLines etc
self.disp(' |1 - info_user/info_numerical| < maxViolation = '+ str(self.maxViolation))
j = -1
for i in range(len(self.input)):
if len(self.input) > 1: self.disp('by input variable number ' + str(i) + ':')
if isinstance(self.d, tuple) and self.d[i] is None:
self.disp('user-provided derivative for input number ' + str(i) + ' is absent, skipping the one;')
self.disp(separator)
continue
if not isinstance(self.input[i], oofun):
self.disp('input number ' + str(i) + ' is not oofun instance, skipping the one;')
self.disp(separator)
continue
j += 1
check_d1(lambda *args: self.fun(*args), ds[j], input, \
func_name=self.name, diffInt=self.diffInt, pointVal = val, args=self.args, \
stencil = max((3, self.stencil)), maxViolation=self.maxViolation, varForCheck = i)
def getOrder(self, Vars=None, fixedVars=None, fixedVarsScheduleID=-1):
# TODO: improve it wrt fixedVarsScheduleID
# returns polinomial order of the oofun
if isinstance(Vars, oofun): Vars = set([Vars])
elif Vars is not None and type(Vars) != set: Vars = set(Vars)
if isinstance(fixedVars, oofun): fixedVars = set([fixedVars])
elif fixedVars is not None and type(fixedVars) != set: fixedVars = set(fixedVars)
sameVarsScheduleID = fixedVarsScheduleID == self._lastOrderVarsID
rebuildFixedCheck = not sameVarsScheduleID
if fixedVarsScheduleID != -1: self._lastOrderVarsID = fixedVarsScheduleID
if rebuildFixedCheck:
# ajust new value of self._order wrt new free/fixed vars schedule
if self.discrete: self._order = 0
if self.is_oovar:
if (fixedVars is not None and self in fixedVars) or (Vars is not None and self not in Vars):
self._order = 0
else:
self._order = 1
else:
# orders = [(inp.getOrder(Vars, fixedVars) if isinstance(inp, oofun) else 0) for inp in self.input]
orders = []
for inp in self.input:
if isinstance(inp, oofun):
orders.append(inp.getOrder(Vars, fixedVars))
elif isinstance(inp, OOArray):
orders += [(elem.getOrder(Vars, fixedVars) if isinstance(elem, oofun) else 0) for elem in inp.view(ndarray)]
self._order = inf if any(asarray(orders) != 0) else 0
return self._order
# TODO: should broadcast return non-void result?
def _broadcast(self, func, useAttachedConstraints, *args, **kwargs):
if self._broadcast_id == oofun._BroadCastID:
return # already done for this one
self._broadcast_id = oofun._BroadCastID
# TODO: possibility of reverse order?
if self.input is not None:
for inp in self.input:
if not isinstance(inp, oofun): continue
inp._broadcast(func, useAttachedConstraints, *args, **kwargs)
if useAttachedConstraints:
for c in self.attachedConstraints:
c._broadcast(func, useAttachedConstraints, *args, **kwargs)
func(self, *args, **kwargs)
# def isUncycled(self):
# # TODO: speedup the function via omitting same oofuns
# # TODO: handle fixed variables
# deps = []
# oofuns = []
# dep_num = []
# for elem in self._getDep():
# if isinstance(elem, oofun):
# if elem.is_oovar:
# deps.append(elem)
# dep_num.append(1)
# else:
# tmp = elem._getDep()
# dep_num.append(len(tmp))
# if tmp is not None:
# deps.append(tmp)
# oofuns.append(elem)
#
# if any([not elem.isUncycled() for elem in oofuns]):
# return False
# r1 = set()
# r1.update(deps)
# return True if len(r1) == sum(dep_num) else False
def uncertainty(self, point, deviations, actionOnAbsentDeviations='warning'):
'''
result = oofun.uncertainty(point, deviations, actionOnAbsentDeviations='warning')
point and deviations should be Python dicts of pairs (oovar, value_for_oovar)
actionOnAbsentDeviations =
'error' (raise FuncDesigner exception) |
'skip' (treat as fixed number with zero deviation) |
'warning' (print warning, treat as fixed number)
Sparse large-scale examples haven't been tested,
we could implement and test it properly on demand
'''
dep = self._getDep()
dev_keys = set(deviations.keys())
set_diff = dep.difference(dev_keys)
nAbsent = len(set_diff)
if actionOnAbsentDeviations != 'skip':
if len(set_diff) != 0:
if actionOnAbsentDeviations == 'warning':
pWarn('''
dict of deviations miss %d variables (oovars): %s;
they will be treated as fixed numbers with zero deviations
''' % (nAbsent, list(set_diff)))
else:
raise FuncDesignerException('dict of deviations miss %d variable(s) (oovars): %s' % (nAbsent, list(set_diff)))
d = self.D(point, exactShape=True) if nAbsent == 0 else self.D(point, fixedVars = set_diff, exactShape=True)
tmp = [dot(val, (deviations[key] if isscalar(deviations[key]) else asarray(deviations[key]).reshape(-1, 1)))**2 for key, val in d.items()]
tmp = [asscalar(elem) if isinstance(elem, ndarray) and elem.size == 1 else elem for elem in tmp]
r = atleast_2d(hstack(tmp)).sum(1)
return r ** 0.5
# For Python 3:
__rtruediv__ = __rdiv__
__truediv__ = __div__
def IMPLICATION(*args, **kw):
raise FuncDesignerException('oofun.IMPLICATION is temporary disabled, use ifThen(...) or IMPLICATION(...) instead')
""" End of class oofun """
# TODO: make it work for ooSystem as well
def broadcast(func, oofuncs, useAttachedConstraints, *args, **kwargs):
if isinstance(oofuncs, oofun):
oofuncs = [oofuncs]
oofun._BroadCastID += 1
for oof in oofuncs:
if oof is not None:
oof._broadcast(func, useAttachedConstraints, *args, **kwargs)
def _getAllAttachedConstraints(oofuns):
from FuncDesigner import broadcast
r = set()
def F(oof):
#print len(oof.attachedConstraints)
r.update(oof.attachedConstraints)
broadcast(F, oofuns, useAttachedConstraints=True)
return r
#def discreteNLH(_input_bool_oofun, Lx, Ux, p, dataType):
#
# T0, res, DefiniteRange = _input_bool_oofun.nlh(Lx, Ux, p, dataType)
# #T = 1.0 - T0
# #R = dict([(v, 1.0-val) for v, val in res.items()])
# return T.flatten(), R, DefiniteRange
def nlh_and(_input, dep, Lx, Ux, p, dataType):
nlh_0 = array(0.0)
R = {}
DefiniteRange = True
elems_nlh = [(elem.nlh(Lx, Ux, p, dataType) if isinstance(elem, oofun) \
else (0, {}, None) if elem is True
else (inf, {}, None) if elem is False
else raise_except()) for elem in _input]
for T0, res, DefiniteRange2 in elems_nlh:
DefiniteRange = logical_and(DefiniteRange, DefiniteRange2)
for T0, res, DefiniteRange2 in elems_nlh:
if T0 is None or T0 is True: continue
if T0 is False or all(T0 == inf):
return inf, {}, DefiniteRange
if all(isnan(T0)):
raise FuncDesignerException('unimplemented for non-oofun or fixed oofun input yet')
if type(T0) == ndarray:
if nlh_0.shape == T0.shape:
nlh_0 += T0
elif nlh_0.size == T0.size:
nlh_0 += T0.reshape(nlh_0.shape)
else:
nlh_0 = nlh_0 + T0
else:
nlh_0 += T0
# debug
# if not any(isfinite(nlh_0)):
# return inf, {}, DefiniteRange
# for T0, res, DefiniteRange2 in elems_nlh:
#debug end
T_0_vect = T0.reshape(-1, 1) if type(T0) == ndarray else T0
for v, val in res.items():
r = R.get(v, None)
if r is None:
R[v] = val - T_0_vect
else:
r += (val if r.shape == val.shape else val.reshape(r.shape)) - T_0_vect
nlh_0_shape = nlh_0.shape
nlh_0 = nlh_0.reshape(-1, 1)
for v, val in R.items():
# TODO: check it
#assert all(isfinite(val))
tmp = val + nlh_0
tmp[isnan(tmp)] = inf # when val = -inf summation with nlh_0 == inf
R[v] = tmp
return nlh_0.reshape(nlh_0_shape), R, DefiniteRange
def nlh_xor(_input, dep, Lx, Ux, p, dataType):
nlh_0 = array(0.0)
nlh_list = []
nlh_list_m = {}
num_inf_m = {}
S_finite = array(0.0)
num_inf_0 = atleast_1d(0)
num_inf_elems = []
R_diff = {}
R_inf = {}
#S_finite_diff = {}
DefiniteRange = True
# elems_nlh = [(elem.nlh(Lx, Ux, p, dataType) if isinstance(elem, oofun) \
# else (0, {}, None) if elem is True
# else (inf, {}, None) if elem is False
# else raise_except()) for elem in _input]
elems_lh = [(elem.lh(Lx, Ux, p, dataType) if isinstance(elem, oofun) \
else (inf, {}, None) if elem is True
else (0, {}, None) if elem is False
else raise_except()) for elem in _input]
for T0, res, DefiniteRange2 in elems_lh:
DefiniteRange = logical_and(DefiniteRange, DefiniteRange2)
for j, (T0, res, DefiniteRange2) in enumerate(elems_lh):
if T0 is None:
raise FuncDesignerException('probably bug in FD kernel')
# !!!!!!!!!!!!!!!!! check "len(elems_lh)" below while calculating P_t
#continue
if all(isnan(T0)):
raise FuncDesignerException('unimplemented for non-oofun or fixed oofun input yet')
#T_0_vect = T0.reshape(-1, 1) if type(T0) == ndarray else T0
T_inf = where(isfinite(T0), 0, 1)
num_inf_elems.append(T_inf)
T0 = where(isfinite(T0), T0, 0.0)
two_pow_t0 = 2.0 ** T0
if type(T0) == ndarray:
if nlh_0.shape == T0.shape:
nlh_0 += T0
num_inf_0 += T_inf
S_finite += two_pow_t0
elif nlh_0.size == T0.size:
nlh_0 += T0.reshape(nlh_0.shape)
num_inf_0 += T_inf.reshape(nlh_0.shape)
S_finite += two_pow_t0.reshape(nlh_0.shape)
else:
nlh_0 = nlh_0 + T0
num_inf_0 = num_inf_0 + T_inf
S_finite = S_finite + two_pow_t0.reshape(nlh_0.shape)
else:
nlh_0 += T0
num_inf_0 += T_inf
S_finite += two_pow_t0
nlh_list.append(T0)
for v, val in res.items():
T_inf_v = where(isfinite(val), 0, 1)
val_noninf = where(isfinite(val), val, 0)
T0v = val_noninf - T0.reshape(-1, 1)
r = nlh_list_m.get(v, None)
if r is None:
nlh_list_m[v] = [(j, T0v)]
num_inf_m[v] = [(j, T_inf_v.copy())]
#num_inf_m[v] = T_inf_v.copy()
else:
r.append((j, T0v))
num_inf_m[v].append((j, T_inf_v.copy()))
#num_inf_m[v] +=T_inf_v
r = R_inf.get(v, None)
T_inf = T_inf.reshape(-1, 1)
if r is None:
R_inf[v] = T_inf_v - T_inf#.reshape(-1, 1)
R_diff[v] = T0v.copy()
else:
# TODO: check for 1st elem of size 1
r += (T_inf_v if r.shape == T_inf_v.shape else T_inf_v.reshape(r.shape)) - T_inf#.reshape(-1, 1)
R_diff[v] += T0v
nlh_1 = [nlh_0 - elem for elem in nlh_list]
# !!! TODO: speedup it via matrix insted of sequence of vectors
num_infs = [num_inf_0 - t for t in num_inf_elems]
S1 = PythonSum([2.0 ** where(num_infs[j] == 0, -t, -inf) for j, t in enumerate(nlh_1)])
S2 = atleast_1d(len(elems_lh) * 2.0 ** (-nlh_0))
S2[num_inf_0 != 0] = 0
#nlh_t = -log(S2 - S1 + 1.0)
#nlh_t = -log1p(S2 - S1) * 1.4426950408889634
nlh_t = -log2(S1-S2)
# assert not any(isnan(nlh_t))
# if not all(isfinite(nlh_t)):
# print('='*10)
# print(nlh_t)
# print(elems_lh)
# raise 0
#print(elems_lh)
# print(R_inf)
#raise 0
R = {}
nlh_0 = nlh_0.reshape(-1, 1)
num_inf_0 = num_inf_0.reshape(-1, 1)
for v, nlh_diff in R_diff.items():
nlh = nlh_0 + nlh_diff
nlh_1 = [nlh - elem.reshape(-1, 1) for elem in nlh_list]
for j, val in nlh_list_m[v]:
nlh_1[j] -= val
Tmp = R_inf[v] + num_inf_0
num_infs = [Tmp] * len(nlh_1)
for j, num_inf in num_inf_m[v]:
num_infs[j] = num_inf
num_infs2 = [Tmp - elem for elem in num_infs]
#num_infs = num_inf - num_inf_m[v]
S1 = PythonSum([2.0 ** where(num_infs2[j] == 0, -elem, -inf) for j, elem in enumerate(nlh_1)])
S2 = atleast_1d(len(elems_lh) * 2.0 ** (-nlh))
S2[Tmp.reshape(S2.shape) != 0] = 0
R[v] = -log2(S1 - S2)
#R[v] = -log1p(S2 - S1) * 1.4426950408889634
for v, val in R.items():
val[isnan(val)] = inf
val[val < 0.0] = 0.0
# print('-'*10)
# #print(Lx, Ux)
# print('elems_lh:', elems_lh)
# print(nlh_t, R, DefiniteRange)
#raw_input()
# if nlh_t.size > 2:
# raise 0
return nlh_t, R, DefiniteRange
#def nlh_not(_input_bool_oofun, dep, Lx, Ux, p, dataType):
# if _input_bool_oofun is True or _input_bool_oofun is False:
# raise 'unimplemented for non-oofun input yet'
# T0, res, DefiniteRange = _input_bool_oofun.nlh(Lx, Ux, p, dataType)
# T = 1.0 - T0
# R = dict([(v, 1.0-val) for v, val in res.items()])
# return T, R, DefiniteRange
def nlh_not(_input_bool_oofun, dep, Lx, Ux, p, dataType):
if _input_bool_oofun is True or _input_bool_oofun is False:
raise 'unimplemented for non-oofun input yet'
T0, res, DefiniteRange = _input_bool_oofun.nlh(Lx, Ux, p, dataType)
T = reverse_l2P(T0)
R = dict([(v, reverse_l2P(val)) for v, val in res.items()])
return T, R, DefiniteRange
def reverse_l2P(l2P):
l2P = atleast_1d(l2P)# elseware bug "0-d arrays cannot be indexed"
#l2P[l2P<0]=0
r = 1.0 / l2P
ind = l2P < 15
r[ind] = -log2(1-2**(-l2P[ind]))
#r[r<0] = 0
return r
def AND(*args):
Args = args[0] if len(args) == 1 and isinstance(args[0], (ndarray, tuple, list, set)) else args
# for arg in Args:
# if isinstance(arg, SmoothFDConstraint) and arg.lb == arg.ub and arg.tol == 0 and not arg.alt_nlh_func:
# pass
#raise FuncDesignerException('equality constraint for smooth func inside logical FD func should have user-assigned tolerance')
assert not isinstance(args[0], ndarray), 'unimplemented yet'
for arg in Args:
if not isinstance(arg, oofun):
raise FuncDesignerException('FuncDesigner logical AND currently is implemented for oofun instances only')
#if other is True: return self
f = logical_and if len(Args) == 2 else alt_AND_engine
r = BooleanOOFun(f, Args, vectorized = True)
r.nlh = lambda *arguments: nlh_and(Args, r._getDep(), *arguments)
r.oofun = r
return r
def alt_AND_engine(*input):
tmp = input[0]
for i in range(1, len(input)):
tmp = logical_and(tmp, input[i])
return tmp
# TODO: multiple args
XOR_prev = lambda arg1, arg2: (arg1 & ~arg2) | (~arg1 & arg2)
def XOR(*args):
Args = args[0] if len(args) == 1 and isinstance(args[0], (ndarray, tuple, list, set)) else args
# for arg in Args:
# if isinstance(arg, SmoothFDConstraint) and arg.lb == arg.ub and arg.tol == 0 and not arg.alt_nlh_func:
# pass
#raise FuncDesignerException('equality constraint for smooth func inside logical FD func should have user-assigned tolerance')
assert not isinstance(args[0], ndarray), 'unimplemented yet'
for arg in Args:
if not isinstance(arg, oofun):
raise FuncDesignerException('FuncDesigner logical XOR currently is implemented for oofun instances only')
#if other is True: return self
#f = lambda *args: logical_xor(hstack([asarray(elem).reshape(-1, 1) for elem in args]))
r = BooleanOOFun(f_xor, Args, vectorized = True)
r.nlh = lambda *arguments: nlh_xor(Args, r._getDep(), *arguments)
r.oofun = r # is it required?
return r
def f_xor(*args):
r = sum(array(args), 0)
return r == 1
EQUIVALENT = lambda arg1, arg2: ((arg1 & arg2) | (~arg1 & ~arg2))
def NOT(_bool_oofun):
assert not isinstance(_bool_oofun, (ndarray, list, tuple, set)), 'disjunctive and other logical constraint are not implemented for ooarrays/ndarrays/lists/tuples yet'
#Args = args[0] if len(args) == 1 and type(args[0]) in (tuple, list, set) else args
#Args = args if type(args) in (tuple, list, set) else [args]
if not isinstance(_bool_oofun, oofun):
raise FuncDesignerException('FuncDesigner logical NOT currently is implemented for oofun instances only')
# if isinstance(_bool_oofun, SmoothFDConstraint) and _bool_oofun.lb == _bool_oofun.ub and _bool_oofun.tol == 0 and not _bool_oofun.alt_nlh_func:
# raise FuncDesignerException('equality constraint for smooth func inside logical FD func should have user-assigned tolerance')
#if other is True: return False
r = BooleanOOFun(logical_not, [_bool_oofun], vectorized = True)
r.oofun = r
if _bool_oofun.is_oovar:
r.lh = lambda *arguments: nlh_not(_bool_oofun, r._getDep(), *arguments)
r.nlh = _bool_oofun.lh
else:
r.nlh = lambda *arguments: nlh_not(_bool_oofun, r._getDep(), *arguments)
return r
NAND = lambda *args, **kw: NOT(AND(*args, **kw))
NOR = lambda *args, **kw: NOT(OR(*args, **kw))
def OR(*args):
Args = args[0] if len(args) == 1 and isinstance(args[0], (ndarray, list, tuple, set)) else args
assert not isinstance(args[0], ndarray), 'unimplemented yet'
for arg in Args:
if not isinstance(arg, oofun):
raise FuncDesignerException('FuncDesigner logical AND currently is implemented for oofun instances only')
r = ~ AND([~elem for elem in Args])
#r.fun = np.logical_or
r.oofun = r
return r
class BooleanOOFun(oofun):
_unnamedBooleanOOFunNumber = 0
discrete = True
#alt_nlh_func = False
# an oofun that returns True/False
def __init__(self, func, _input, *args, **kwargs):
oofun.__init__(self, func, _input, *args, **kwargs)
#self.input = oofun_Involved.input
BooleanOOFun._unnamedBooleanOOFunNumber += 1
self.name = 'unnamed_boolean_oofun_' + str(BooleanOOFun._unnamedBooleanOOFunNumber)
self.oofun = oofun(lambda *args, **kw: asanyarray(func(*args, **kw), int8), _input, vectorized = True)
# TODO: THIS SHOULD BE USED IN UP-LEVEL ONLY
self.lb = self.ub = 1
__hash__ = lambda self: self._id
def size(self, *args, **kwargs): raise FuncDesignerException('currently BooleanOOFun.size() is disabled')
def D(self, *args, **kwargs): raise FuncDesignerException('currently BooleanOOFun.D() is disabled')
def _D(self, *args, **kwargs): raise FuncDesignerException('currently BooleanOOFun._D() is disabled')
def nlh(self, *args, **kw):
raise FuncDesignerException('This is virtual method to be overloaded in derived class instance')
__and__ = AND
#IMPLICATION = IMPLICATION
__eq__ = EQUIVALENT
__ne__ = lambda self, arg: NOT(self==arg)
def __or__(self, other):
#if other is False: return self
return ~((~self) & (~other))
# print('__or__')
# r = BooleanOOFun(logical_or, (self, other), vectorized = True)
# r.nlh = lambda *args: nlh_or((self, other), r._getDep(), *args)
# r.oofun = r
# return r
def __xor__(self, other):
return BooleanOOFun(logical_xor, (self, other), vectorized = True)
def __invert__(self):
r = BooleanOOFun(logical_not, self, vectorized = True)
r.nlh = lambda *args: nlh_not(self, r._getDep(), *args)
r.oofun = r
return r
class BaseFDConstraint(BooleanOOFun):
isConstraint = True
tol = 0.0
expected_kwargs = set(['tol', 'name'])
#def __getitem__(self, point):
def __call__(self, *args, **kwargs):
expected_kwargs = self.expected_kwargs
if not set(kwargs.keys()).issubset(expected_kwargs):
raise FuncDesignerException('Unexpected kwargs: should be in '+str(expected_kwargs)+' got: '+str(kwargs.keys()))
for elem in expected_kwargs:
if elem in kwargs:
setattr(self, elem, kwargs[elem])
if len(args) > 1: raise FuncDesignerException('No more than single argument is expected')
if len(args) == 0:
if len(kwargs) == 0: raise FuncDesignerException('You should provide at least one argument')
return self
if isinstance(args[0], str):
self.name = args[0]
return self
elif hasattr(args[0], 'xf'):
return self(args[0].xf)
return self._getFuncCalcEngine(*args, **kwargs)
def _getFuncCalcEngine(self, *args, **kwargs):
if not isinstance(args[0], dict):
raise FuncDesignerException('unexpected type: %s' % type(args[0]))
isMultiPoint = isinstance(args[0], ooPoint) and args[0].isMultiPoint == True
val = self.oofun(args[0])
Tol = max((0.0, self.tol))
if isMultiPoint:
return logical_and(self.lb-Tol<= val, val <= self.ub + Tol)
elif any(isnan(val)):
return False
if any(atleast_1d(self.lb-val)>Tol):
return False
elif any(atleast_1d(val-self.ub)>Tol):
return False
return True
def __init__(self, oofun_Involved, *args, **kwargs):
BooleanOOFun.__init__(self, oofun_Involved._getFuncCalcEngine, (oofun_Involved.input if not oofun_Involved.is_oovar else oofun_Involved), *args, **kwargs)
#oofun.__init__(self, lambda x: oofun_Involved(x), input = oofun_Involved)
if len(args) != 0:
raise FuncDesignerException('No args are allowed for FuncDesigner constraint constructor, only some kwargs')
# TODO: replace self.oofun by self.engine
self.oofun = oofun_Involved
class SmoothFDConstraint(BaseFDConstraint):
__getitem__ = lambda self, point: self.__call__(point)
def __init__(self, *args, **kwargs):
BaseFDConstraint.__init__(self, *args, **kwargs)
self.lb, self.ub = -inf, inf
for key, val in kwargs.items():
if key in ['lb', 'ub', 'tol']:
setattr(self, key, asfarray(val))
else:
raise FuncDesignerException('Unexpected key in FuncDesigner constraint constructor kwargs')
def lh(self, *args, **kw): # overwritten in ooVar, mb something else
if '_invert' not in self.__dict__:
self._invert = NOT(self)
return self._invert.nlh(*args, **kw)
def nlh(self, Lx, Ux, p, dataType):
m = Lx.shape[0] # is changed in the cycle
if m == 0:
assert 0, 'bug in FuncDesigner'
#return None, None, None
DefiniteRange = True
tol = self.tol if self.tol > 0.0 else p.contol if self.tol == 0 else 0.0 # 0 for negative tolerances
# TODO: check it
if p.solver.dataHandling == 'sorted': tol = 0
selfDep = (self.oofun._getDep() if not self.oofun.is_oovar else set([self.oofun]))
# prev
#domainData = [(v, (Lx[:, k], Ux[:, k])) for k, v in enumerate(p._freeVarsList)]
# new
# TODO: improve it
domainData = [(v, (Lx[:, k], Ux[:, k])) for k, v in enumerate(p._freeVarsList) if v in selfDep]
domain = ooPoint(domainData, skipArrayCast=True)
domain.isMultiPoint = True
domain.dictOfFixedFuncs = p.dictOfFixedFuncs
r, r0 = self.oofun.iqg(domain, dataType, self.lb, self.ub)
Lf, Uf = r0.lb, r0.ub
tmp = getSmoothNLH(tile(Lf, (2, 1)), tile(Uf, (2, 1)), self.lb, self.ub, tol, m, dataType)
T02 = tmp
T0 = T02[:, tmp.shape[1]/2:].flatten()
res = {}
if len(r):
dep = selfDep.intersection(domain.keys()) # TODO: Improve it
for v in dep:
Lf, Uf = vstack((r[v][0].lb, r[v][1].lb)), vstack((r[v][0].ub, r[v][1].ub))
# TODO: 1) FIX IT it for matrix definiteRange
# 2) seems like DefiniteRange = (True, True) for any variable is enough for whole range to be defined in the involved node
DefiniteRange = logical_and(DefiniteRange, r[v][0].definiteRange)
DefiniteRange = logical_and(DefiniteRange, r[v][1].definiteRange)
tmp = getSmoothNLH(Lf, Uf, self.lb, self.ub, tol, m, dataType) #- T02
#tmp[isnan(tmp)] = inf
res[v] = tmp
return T0, res, DefiniteRange
def getSmoothNLH(Lf, Uf, lb, ub, tol, m, dataType):
M = prod(Lf.shape) / (2*m)
#init
Lf, Uf = Lf.reshape(2*M, m).T, Uf.reshape(2*M, m).T
#1
#Lf, Uf = Lf.reshape(m, 2*M), Uf.reshape(m, 2*M)
lf1, lf2, uf1, uf2 = Lf[:, 0:M], Lf[:, M:2*M], Uf[:, 0:M], Uf[:, M:2*M]
UfLfDiff = Uf - Lf
if UfLfDiff.dtype.type in [int8, int16, int32, int64, int]:
UfLfDiff = asfarray(UfLfDiff)
#UfLfDiff[UfLfDiff > 1e200] = 1e200
if lb == ub:
val = ub
ind1, ind2 = val - tol > Uf, val+tol < Lf
# residual[ind1] += val - tol - Uf[ind1]
# residual[ind2] += Lf[ind2] - (val + tol)
Uf_t, Lf_t = Uf.copy(), Lf.copy()
if Uf.dtype.type in [int8, int16, int32, int64, int] or Lf.dtype.type in [int8, int16, int32, int64, int]:
Uf_t, Lf_t = asfarray(Uf_t), asfarray(Lf_t)
Uf_t[Uf_t > val + tol] = val + tol
Lf_t[Lf_t < val - tol] = val - tol
allowedLineSegmentLength = Uf_t - Lf_t
tmp = allowedLineSegmentLength / UfLfDiff
#tmp[tmp>1] = 1
tmp[logical_or(isinf(Lf), isinf(Uf))] = 1e-10 # (to prevent inf/inf=nan); TODO: rework it
tmp[allowedLineSegmentLength == 0.0] = 1.0 # may be encountered if Uf == Lf, especially for integer probs
tmp[tmp<1e-300] = 1e-300 # TODO: improve it
# TODO: for non-exact interval quality increase nlh while moving from 0.5*(Ux-Lx)
tmp[val - tol > Uf] = 0
tmp[val + tol < Lf] = 0
elif isfinite(lb) and not isfinite(ub):
tmp = (Uf - (lb - tol)) / UfLfDiff
#ind = Lf < lb-tol
#residual[ind] += lb-Lf[ind]
tmp[logical_and(isinf(Lf), logical_not(isinf(Uf)))] = 1e-10 # (to prevent inf/inf=nan); TODO: rework it
tmp[isinf(Uf)] = 1-1e-10 # (to prevent inf/inf=nan); TODO: rework it
tmp[tmp<1e-300] = 1e-300 # TODO: improve it
tmp[tmp>1.0] = 1.0
#tmp[lb-tol> Uf] = 0
tmp[lb-tol > Uf] = 0
tmp[lb <= Lf] = 1
#tmp[lb <= Lf] = 1
elif isfinite(ub) and not isfinite(lb):
tmp = (ub + tol - Lf ) / UfLfDiff
#ind = Uf > ub+tol
#residual[ind] += Uf[ind]-ub
tmp[isinf(Lf)] = 1-1e-10 # (to prevent inf/inf=nan);TODO: rework it
tmp[logical_and(isinf(Uf), logical_not(isinf(Lf)))] = 1e-10 # (to prevent inf/inf=nan); TODO: rework it
tmp[tmp<1e-300] = 1e-300 # TODO: improve it
tmp[tmp>1.0] = 1.0
tmp[ub+tol < Lf] = 0
#tmp[ub < Lf] = 0
tmp[ub >= Uf] = 1
#tmp[ub >= Uf] = 1
else:
raise FuncDesignerException('this part of interalg code is unimplemented for double-box-bound constraints yet')
tmp = -log2(tmp)
tmp[isnan(tmp)] = inf # to prevent some issues in disjunctive cons
return tmp
class Constraint(SmoothFDConstraint):
def __init__(self, *args, **kwargs):
SmoothFDConstraint.__init__(self, *args, **kwargs)
class BoxBoundConstraint(SmoothFDConstraint):
def __init__(self, *args, **kwargs):
SmoothFDConstraint.__init__(self, *args, **kwargs)
class Derivative(dict):
def __init__(self):
pass
#def ooFun(*args, **kwargs):
# r = oofun(*args, **kwargs)
# r.isCostly = True
# return r
def atleast_oofun(arg):
if isinstance(arg, oofun):
return arg
elif hasattr(arg, 'copy'):
tmp = arg.copy()
return oofun(lambda *args, **kwargs: tmp, input = None, getOrder = lambda *args, **kwargs: 0, discrete=True)#, isConstraint = True)
elif isscalar(arg):
tmp = array(arg, 'float')
return oofun(lambda *args, **kwargs: tmp, input = None, getOrder = lambda *args, **kwargs: 0, discrete=True)#, isConstraint = True)
else:
#return oofun(lambda *args, **kwargs: arg(*args, **kwargs), input=None, discrete=True)
raise FuncDesignerException('incorrect type for the function _atleast_oofun')
def mul_aux_d(x, y):
Xsize, Ysize = Len(x), Len(y)
if Xsize == 1:
return Copy(y)
elif Ysize == 1:
return Diag(None, scalarMultiplier = y, size = Xsize)
elif Xsize == Ysize:
return Diag(y)
else:
raise FuncDesignerException('for oofun multiplication a*b should be size(a)=size(b) or size(a)=1 or size(b)=1')
def mul_interval(self, other, isOtherOOFun, domain, dtype):#*args, **kw):
if domain.isMultiPoint and isOtherOOFun and self.is_oovar and (self.domain is bool or self.domain is 'bool'):
#ind_nz = where(domain[self][1]!=0)[0]
lb_ub, definiteRange = other._interval(domain, dtype)
n = domain[self][1].size
R = zeros((2, n), dtype)
ind = domain[self][0]!=0
R[0][ind] = lb_ub[0][ind]
ind = domain[self][1]!=0
R[1][ind] = lb_ub[1][ind]
return R, definiteRange
# if ind_nz.size == 0:
# n = domain[self][1].size
# R = zeros((2, n), dtype)
# definiteRange = empty(n, bool)
# definiteRange.fill(True)# TODO: sometimes it may be wrong
# return R, definiteRange
# TODO: implement it
# subdomain = ooPoint([(v, (d[0][ind_nz], d[1][ind_nz])) for v, d in domain.items()], skipArrayCast=True)
# subdomain.isMultiPoint = True
# #subdomain.modificationVar = domain.modificationVar
# #subdomain.useSave = domain.useSave
# #subdomain.localStoredIntervals = {}
# lb_ub, definiteRange2 = other._interval(subdomain, dtype)
# R[:, ind_nz] = lb_ub.copy()
# definiteRange[ind_nz] = definiteRange2
# return R, definiteRange
#
# elif isOtherOOFun:
# # TODO: implement it
# pass
#
# # changes end
lb1_ub1, definiteRange = self._interval(domain, dtype)
lb1, ub1 = lb1_ub1[0], lb1_ub1[1]
if isOtherOOFun:
lb2_ub2, definiteRange2 = other._interval(domain, dtype)
definiteRange = logical_and(definiteRange, definiteRange2)
lb2, ub2 = lb2_ub2[0], lb2_ub2[1]
else:
lb2, ub2 = other, other # TODO: improve it
if all(lb2 >= 0) and all(lb1 >= 0):
t_min, t_max = atleast_1d(lb1 * lb2), atleast_1d(ub1 * ub2)
elif isscalar(other):
t_min, t_max = (lb1 * other, ub1 * other) if other >= 0 else (ub1 * other, lb1 * other)
else:
if isOtherOOFun:
t = vstack((lb1 * lb2, ub1 * lb2, \
lb1 * ub2, ub1 * ub2))# TODO: improve it
else:
t = vstack((lb1 * other, ub1 * other))# TODO: improve it
t_min, t_max = atleast_1d(nanmin(t, 0)), atleast_1d(nanmax(t, 0))
#assert isinstance(t_min, ndarray) and isinstance(t_max, ndarray), 'Please update numpy to more recent version'
if any(isinf(lb1)) or any(isinf(lb2)) or any(isinf(ub1)) or any(isinf(ub2)):
ind1_zero_minus = logical_and(lb1<0, ub1>=0)
ind1_zero_plus = logical_and(lb1<=0, ub1>0)
ind2_zero_minus = logical_and(lb2<0, ub2>=0)
ind2_zero_plus = logical_and(lb2<=0, ub2>0)
has_plus_inf_1 = logical_or(logical_and(ind1_zero_minus, lb2==-inf), logical_and(ind1_zero_plus, ub2==inf))
has_plus_inf_2 = logical_or(logical_and(ind2_zero_minus, lb1==-inf), logical_and(ind2_zero_plus, ub1==inf))
# !!!! lines with zero should be before lines with inf !!!!
ind = logical_or(logical_and(lb1==-inf, ub2==0), logical_and(lb2==-inf, ub1==0))
t_max[atleast_1d(logical_and(ind, t_max<0.0))] = 0.0
t_max[atleast_1d(logical_or(has_plus_inf_1, has_plus_inf_2))] = inf
t_max[atleast_1d(logical_or(logical_and(lb1==0, ub2==inf), logical_and(lb2==0, ub1==inf)))] = inf
has_minus_inf_1 = logical_or(logical_and(ind1_zero_plus, lb2==-inf), logical_and(ind1_zero_minus, ub2==inf))
has_minus_inf_2 = logical_or(logical_and(ind2_zero_plus, lb1==-inf), logical_and(ind2_zero_minus, ub1==inf))
# !!!! lines with zero should be before lines with -inf !!!!
t_min[atleast_1d(logical_or(logical_and(lb1==0, ub2==inf), logical_and(lb2==0, ub1==inf)))] = 0.0
t_min[atleast_1d(logical_or(logical_and(lb1==-inf, ub2==0), logical_and(lb2==-inf, ub1==0)))] = 0.0
t_min[atleast_1d(logical_or(has_minus_inf_1, has_minus_inf_2))] = -inf
# assert not any(isnan(t_min)) and not any(isnan(t_max))
return vstack((t_min, t_max)), definiteRange
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/ooPoint.py 000644 000766 000024 00000005710 12100431107 025411 0 ustar 00rfadams staff 000000 000000 # created by Dmitrey
#from numpy import inf, asfarray, copy, all, any, empty, atleast_2d, zeros, dot, asarray, atleast_1d, empty, ones, ndarray, \
#where, array, nan, ix_, vstack, eye, array_equal, isscalar, diag, log, hstack, sum, prod, nonzero, isnan
#from numpy.linalg import norm
#from misc import FuncDesignerException, Diag, Eye, pWarn, scipyAbsentMsg
#from copy import deepcopy
from FDmisc import FuncDesignerException
from baseClasses import Stochastic
from numpy import asanyarray, ndarray, isscalar
try:
from scipy.sparse import isspmatrix
except ImportError:
isspmatrix = lambda *args, **kw: False
Len = lambda x: 1 if isscalar(x) else x.size if type(x)==ndarray else len(x)
def ooMultiPoint(*args, **kw):
kw['skipArrayCast'] = True
r = ooPoint(*args, **kw)
r.isMultiPoint = True
return r
class ooPoint(dict):
_id = 0
isMultiPoint = False
modificationVar = None # default: no modification variable
useSave = True
useAsMutable = False
def __init__(self, *args, **kwargs):
self.storedIntervals = {}
self.storedSums = {}
self.dictOfFixedFuncs = {}
for fn in ('isMultiPoint', 'modificationVar', 'useSave', 'useAsMutable', 'maxDistributionSize'):
tmp = kwargs.get(fn, None)
if tmp is not None:
setattr(self, fn, tmp)
if kwargs.get('skipArrayCast', False):
Asanyarray = lambda arg: arg
else:
Asanyarray = lambda arg: asanyarray(arg) if not isinstance(arg, Stochastic) else arg#if not isspmatrix(arg) else arg
# TODO: remove float() after Python 3 migraion
if args:
if not isinstance(args[0], dict):
items = [(key, Asanyarray(val) if not isscalar(val) else float(val) if type(val) == int else val) for key, val in args[0]]
else:
items = [(key, Asanyarray(val) if not isscalar(val) else float(val) if type(val) == int else val) for key, val in args[0].items()]
elif kwargs:
items = [(key, Asanyarray(val) if not isscalar(val) else float(val) if type(val) == int else val) for key, val in kwargs.items()]
else:
raise FuncDesignerException('incorrect oopoint constructor arguments')
dict.__init__(self, items)
# TODO: fix it wrt ode2.py
# for key, val in items:
# #assert type(val) not in [list, ndarray] or type(val[0]) != int
# if 'size' in key.__dict__ and type(key.size) == int and Len(val) != key.size:
# s = 'incorrect size for oovar %s: %d is required, %d is obtained' % (key, self.size, Size)
# raise FuncDesignerException(s)
ooPoint._id += 1
self._id = ooPoint._id
def __setitem__(self, *args, **kwargs):
if not self.useAsMutable:
raise FuncDesignerException('ooPoint must be immutable')
dict.__setitem__(self, *args, **kwargs)
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/ooSystem.py 000644 000766 000024 00000023646 12100431107 025614 0 ustar 00rfadams staff 000000 000000 from FDmisc import FuncDesignerException, pWarn
from ooFun import oofun, BaseFDConstraint, _getAllAttachedConstraints
from ooarray import ooarray
from ooPoint import ooPoint
from numpy import isnan, ndarray, isfinite, asscalar, all, asarray, atleast_1d, array_equal
from sle import sle
class ooSystem:
def __init__(self, *args, **kwargs):
assert len(kwargs) == 0, 'ooSystem constructor has no implemented kwargs yet'
self.items = set()
self.constraints = set()
self.nlpSolvers = ['ralg']
self.lpSolvers = ['glpk', 'lpSolve', 'cvxopt_lp', 'cplex','pclp','nlp:ipopt', 'nlp:algencan', 'nlp:scipy_slsqp', 'nlp:ralg']
for arg in args:
if isinstance(arg, set):
self.items.update(arg)
elif isinstance(arg, list) or isinstance(arg, tuple):
self.items.update(set(arg))
elif isinstance(arg, oofun):
self.items.add(arg)
else:
raise FuncDesignerException('Incorrect type %s in ooSystem constructor' % type(arg))
self._changed = True
# [] yields result wrt allowed contol:
#__getitem__ = lambda self, point: ooSystemState([(elem, elem[point]) for elem in self])
def __getitem__(self, item):
raise FuncDesignerException('ooSystem __getitem__ is reserved for future purposes')
def __iadd__(self, *args, **kwargs):
assert len(kwargs) == 0, 'not implemented yet'
self._changed = True
if type(args[0]) in [list, tuple, set]:
assert len(args) == 1
Args = args[0]
else:
Args = args
for elem in Args:
if not isinstance(elem, oofun): raise FuncDesignerException('ooSystem operation += expects only oofuns')
self.items.update(set(Args))
return self
def __iand__(self, *args, **kwargs):
assert len(kwargs) == 0, 'not implemented yet'
self._changed = True
if type(args[0]) in [list, tuple, set]:
assert len(args) == 1
Args = args[0]
else:
Args = args
for elem in Args:
if isinstance(elem, (ooarray, list, tuple)):
for tmp_elem in elem: self.__iand__(tmp_elem)
else:
if not isinstance(elem, BaseFDConstraint): raise FuncDesignerException('ooSystem operation &= expects only FuncDesigner constraints')
#self.constraints.update(set(Args))
self.constraints.add(elem)
return self
# TODO: provide a possibility to yield exact result (ignoring contol)
def __call__(self, point):
if not isinstance(point, dict): raise FuncDesignerException('argument should be Python dictionary')
if not isinstance(point, ooPoint):
point = ooPoint(point)
r = ooSystemState([(elem, simplify(elem(point))) for elem in self.items])
# handling constraints
#!!!!!!!!!!!!!!!!!!! TODO: perform attached constraints lookup only once if ooSystem wasn't modified by += or &= etc
cons = self._getAllConstraints()
activeConstraints = []
allAreFinite = all([all(isfinite(asarray(elem(point)))) for elem in self.items])
for c in cons:
val = c.oofun(point)
if c(point) is False or any(isnan(atleast_1d(val))):
activeConstraints.append(c)
#_activeConstraints.append([c, val, max((val-c.ub, c.lb-val)), c.tol])
r.isFeasible = True if len(activeConstraints) == 0 and allAreFinite else False
#r._activeConstraints = activeConstraints
r.activeConstraints = activeConstraints
return r
def minimize(self, *args, **kwargs):
kwargs['goal'] = 'minimum'
return self._solve(*args, **kwargs)
def maximize(self, *args, **kwargs):
kwargs['goal'] = 'maximum'
return self._solve(*args, **kwargs)
def solve(self, *args, **kwargs):
kwargs['goal'] = 'solution'
assert len(args) == 0 or all([not isinstance(arg, oofun) for arg in args])
return self._solve(*args, **kwargs)
def _solve(self, *args, **kwargs):
#!!!!! TODO: determing LP, MILP, MINLP if possible
try:
import openopt
except ImportError:
raise FuncDesignerException('to perform the operation you should have openopt installed')
constraints = self._getAllConstraints()
# TODO: mention it in doc
if 'constraints' in kwargs:
tmp = set(kwargs['constraints'])
tmp.update(set(constraints))
kwargs['constraints'] = tmp
else:
kwargs['constraints'] = constraints
freeVars, fixedVars = kwargs.get('freeVars', None), kwargs.get('fixedVars', None)
isSystemOfEquations = kwargs['goal'] == 'solution'
isLinear = all([(c.oofun.getOrder(freeVars, fixedVars) < 2) for c in constraints])
if isSystemOfEquations:
if isLinear: # Linear equations system
p = sle(list(kwargs['constraints']), *args, **kwargs)
else: # Nonlinear equations system
f = kwargs['constraints']
############################################
#cons = self._getAllConstraints()
# OLD
# if not all([array_equal(c.lb, c.ub) for c in f]):
# raise FuncDesignerException('Solving constrained nonlinear systems is not implemented for oosystem yet, use SNLE constructor instead')
# kwargs['constraints'] = ()
# NEW
C, F = [], []
for c in f:
F.append(c) if array_equal(c.lb, c.ub) else C.append(c)
kwargs['constraints'] = C
f = F
############################################
p = openopt.SNLE(f, *args, **kwargs)
if 'nlpSolver' in kwargs:
p.solver = kwargs['nlpSolver']
else: # an optimization problem
assert len(args) > 0, 'you should have objective function as 1st argument'
objective = args[0]
if isinstance(objective, BaseFDConstraint):
raise FuncDesignerException("1st argument can't be of type 'FuncDesigner constraint', it should be 'FuncDesigner oofun'")
elif not isinstance(objective, oofun):
raise FuncDesignerException('1st argument should be objective function of type "FuncDesigner oofun"')
isLinear &= objective.getOrder(freeVars, fixedVars) < 2
if isLinear:
p = openopt.LP(*args, **kwargs)
if 'solver' not in kwargs:
for solver in self.lpSolvers:
if (':' not in solver and not openopt.oosolver(solver).isInstalled )or (solver == 'glpk' and not openopt.oosolver('cvxopt_lp').isInstalled):
continue
if solver == 'glpk' :
p2 = openopt.LP([1, -1], lb = [1, 1], ub=[10, 10])
try:
r = p2.solve('glpk', iprint=-5)
except:
continue
if r.istop < 0:
continue
else:
break
if ':' in solver:
pWarn('You have linear problem but no linear solver (lpSolve, glpk, cvxopt_lp) is installed; converter to NLP will be used.')
p.solver = solver
else:
solverName = kwargs['solver']
if type(solverName) != str:
solverName = solverName.__name__
if solverName not in self.lpSolvers:
solverName = 'nlp:' + solverName
p.solver = solverName
p.warn('you are solving linear problem with non-linear solver')
else:
p = openopt.NLP(*args, **kwargs)
if 'solver' not in kwargs:
p.solver = 'ralg'
# TODO: solver autoselect
#if p.iprint >= 0: p.disp('The optimization problem is ' + p.probType)
p._isFDmodel = lambda *args, **kwargs: True
return p.solve() if kwargs.get('manage', False) in (False, 0) else p.manage()
def _getAllConstraints(self):
if self._changed:
cons = self.constraints
cons.update(_getAllAttachedConstraints(self.items | self.constraints))
self._AllConstraints = cons
self._changed = False
return self._AllConstraints
####################### ooSystemState #########################
class ooSystemState:
def __init__(self, keysAndValues, *args, **kwargs):
assert len(args) ==0
assert len(kwargs) ==0
#dict.__init__(self, *args, **kwargs)
self._byID = dict([(key, val) for key, val in keysAndValues])
self._byNames = dict([(key.name, val) for key, val in keysAndValues])
#self.update(self._byNames)
__repr__ = lambda self: ''.join(['\n'+key+' = '+str(val) for key, val in self._byNames.items()])[1:]
def __call__(self, *args, **kwargs):
assert len(kwargs) == 0, "ooSystemState method '__call__' has no implemented kwargs yet"
r = [(self._byNames[arg] if isinstance(arg, str) else self._byID[arg]) for arg in args]
return r[0] if len(r)==1 else r
#__call__ = _call
#__getitem__ = _call
# TODO: implement it
# access by elements and their names
simplify = lambda val: asscalar(val) if isinstance(val, ndarray) and val.size == 1 else val
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/ooVar.py 000644 000766 000024 00000035601 12100431107 025052 0 ustar 00rfadams staff 000000 000000 # created by Dmitrey
from numpy import asarray, empty, inf, any, array, \
asfarray, isscalar, ndarray, int16, int32, int64, float64, tile, vstack, searchsorted, logical_or, where, \
asanyarray, arange, log2, logical_and, ceil
import numpy as np
from FDmisc import FuncDesignerException, isPyPy
from ooFun import oofun, BooleanOOFun, AND, OR, NOT, EQUIVALENT
#from FuncDesigner import IMPLICATION
from ooarray import ooarray
#from FuncDesigner.Interval import adjust_lx_WithDiscreteDomain, adjust_ux_WithDiscreteDomain
from baseClasses import Stochastic
f_none = lambda *args, **kw: None
class oovar(oofun):
is_oovar = True
domain = None
lb = -inf
ub = inf
#shape = nan
#fixed = False
#initialized = False
_unnamedVarNumber = 1#static variable for oovar class
__hash__ = oofun.__hash__
def __init__(self, name=None, *args, **kwargs):
if len(args) > 0: raise FuncDesignerException('incorrect args number for oovar constructor')
if name is None:
self.name = 'unnamed_' + str(oovar._unnamedVarNumber)
oovar._unnamedVarNumber += 1
else:
kwargs['name'] = name
oofun.__init__(self, f_none, *args, **kwargs)
def _interval_(self, domain, dtype = float64):
tmp = domain.get(self, None)
if tmp is None: return None
if isinstance(tmp, ndarray) or isscalar(tmp): # thus variable value is fixed for this calculation
tmp = asarray(tmp, dtype)
return tile(tmp, (2, 1)), True
infinum, supremum = tmp
if type(infinum) in (list, tuple):
infinum = array(infinum, dtype)
elif isscalar(infinum):
infinum = dtype(infinum)
if type(supremum) in (list, tuple):
supremum = array(supremum, dtype)
elif isscalar(supremum):
supremum = dtype(supremum)
# if modificationVar is self:
# assert dtype in (float, float64), 'other types unimplemented yet'
# middle = 0.5 * (supremum+infinum)
# return (vstack((infinum, middle)), True), (vstack((middle, supremum)), True)
return vstack((infinum, supremum)), True
#def _interval_(self, domain, dtype):
#return self._interval()
#_interval_ = _interval
def _getFuncCalcEngine(self, x, **kwargs):
if hasattr(x, 'xf'):
#return x.xf[self]
if x.probType == 'MOP':
s = 'evaluation of MOP result on arguments is unimplemented yet, use r.solutions'
raise FuncDesignerException(s)
return self._getFuncCalcEngine(x.xf, **kwargs) # essential for SP
r = x.get(self, None)
if r is not None:
if isinstance(r, Stochastic):
sz = getattr(x, 'maxDistributionSize', 0)
if sz == 0:
s = '''
if one of function arguments is stochastic distribution
without resolving into quantified value
(e.g. uniform(-10,10) instead of uniform(-10,10, 100), 100 is number of point to emulate)
then you should evaluate the function
onto oopoint with assigned parameter maxDistributionSize'''
raise FuncDesignerException(s)
if not r.quantified:
r = r._yield_quantified(sz)
r = r.copy()
r.stochDep = {self:1}
r.maxDistributionSize = sz
if r.size > sz:
r.reduce(sz)
tmp = getattr(x, '_p', None)
if tmp is not None:
r._p = tmp
return r
r = x.get(self.name, None)
if r is not None:
return r
else:
s = '''for oovar %s the point involved doesn't contain
neither name nor the oovar instance.
Maybe you try to get function value or derivative
in a point where value for an oovar is missing
or run optimization problem
without setting initial value for this variable in start point
''' % self.name
raise FuncDesignerException(s)
def nlh(self, Lx, Ux, p, dataType, other=None):
T0, res, DefiniteRange = get_P(self, Lx, Ux, p, dataType, other, goal_is_nlh = True)
if type(T0) == bool:
assert len(res) == 0
return T0, {}, DefiniteRange
else:
return T0, {self: res}, DefiniteRange
def lh(self, Lx, Ux, p, dataType, other=None):
#print('lh')
T0, res, DefiniteRange = get_P(self, Lx, Ux, p, dataType, other, goal_is_nlh = False)
if type(T0) == bool:
assert len(res) == 0
return T0, {}, DefiniteRange
else:
return T0, {self: res}, DefiniteRange
__and__ = AND
__or__ = OR
#implication = IMPLICATION
__invert__ = NOT
__ne__ = lambda self, arg: NOT(self==arg)
def __eq__(self, other):
if (self.domain is bool or self.domain is 'bool') and isinstance(other, (oovar, BooleanOOFun)):
return EQUIVALENT(self, other)
else:
return oofun.__eq__(self, other)
def formAuxDomain(self):
if 'aux_domain' in self.__dict__: return
self.domain = asanyarray(self.domain)
d = self.domain
# if d.dtype.type not in [string_, unicode, str]:
# raise FuncDesignerException('to compare string with oovar latter should have domain of string type')
if any(d[1:] < d[:-1]):
d.sort()
#self.ub = d.size - 1
D = int(2 ** ceil(log2(d.size)))
self.domain, self.aux_domain = arange(D), d
# if isinstance(x, dict):
# tmp = x.get(self, None)
# if tmp is not None:
# return tmp #if type(tmp)==ndarray else asfarray(tmp)
# elif self.name in x:
# return asfarray(x[self.name])
# else:
# s = 'for oovar ' + self.name + \
# " the point involved doesn't contain niether name nor the oovar instance. Maybe you try to get function value or derivative in a point where value for an oovar is missing"
# raise FuncDesignerException(s)
# elif hasattr(x, 'xf'):
# # TODO: possibility of squeezing
# return x.xf[self]
# else:
# raise FuncDesignerException('Incorrect data type (%s) while obtaining oovar %s value' %(type(x), self.name))
# def _initialize(self, p):
#
# """ Handling size and shape """
# sizes = set()
# shapes = set()
# for fn in ['v0', 'lb', 'ub']:
# if hasattr(self, fn):
# setattr(self, fn, asarray(getattr(self, fn)))
# shapes.add(getattr(self, fn).shape)
# sizes.add(getattr(self, fn).size)
# if self.shape is not nan:
# shapes.add(self.shape)
# sizes.add(prod(self.shape))
# if self.size is not nan: sizes.add(self.size)
# #if len(shapes) > 1: p.err('for oovar fields (if present) lb, ub, v0 should have same shape')
# #elif len(shapes) == 1: self.shape = shapes.pop()
# if len(shapes) >= 1: self.shape = prod(shapes.pop())
#
# if len(sizes) > 1: p.err('for oovar fields (if present) lb, ub, v0 should have same size')
# elif len(sizes)==1 : self.size = sizes.pop()
#
# if self.size is nan: self.size = asarray(self.shape).prod()
# if self.shape is nan:
# assert isfinite(self.size)
# self.shape = (self.size, )
#
#
# """ Handling init value """
## if not hasattr(self, 'lb'):
## self.lb = empty(self.shape)
## self.lb.fill(-inf)
## if not hasattr(self, 'ub'):
## self.ub = empty(self.shape)
## self.ub.fill(inf)
## if any(self.lb > self.ub):
## p.err('lower bound exceeds upper bound, solving impossible')
# if not hasattr(self, 'v0'):
# #p.warn('got oovar w/o init value')
# v0 = zeros(self.shape)
#
# ind = isfinite(self.lb) & isfinite(self.ub)
# v0[ind] = 0.5*(self.lb[ind] + self.ub[ind])
#
# ind = isfinite(self.lb) & ~isfinite(self.ub)
# v0[ind] = self.lb[ind]
#
# ind = ~isfinite(self.lb) & isfinite(self.ub)
# v0[ind] = self.ub[ind]
#
# self.v0 = v0
#
# self.initialized = True
def oovars(*args, **kw):
if isPyPy:
raise FuncDesignerException('''
for PyPy using oovars() is impossible yet.
You could use oovar(size=n), also
you can create list or tuple of oovars in a cycle, e.g.
a = [oovar('a'+str(i)) for i in range(100)]
but you should ensure you haven't operations like k*a or a+val in your code,
it may work in completely different way (e.g. k*a will produce Python list of k a instances)
''')
lb = kw.pop('lb', None)
ub = kw.pop('ub', None)
if len(args) == 1:
if type(args[0]) in (int, int16, int32, int64):
r = ooarray([oovar(**kw) for i in range(args[0])])
elif type(args[0]) in [list, tuple]:
r = ooarray([oovar(name=args[0][i], **kw) for i in range(len(args[0]))])
elif type(args[0]) == str:
r = ooarray([oovar(name=s, **kw) for s in args[0].split()])
else:
raise FuncDesignerException('incorrect args number for oovars constructor')
else:
r = ooarray([oovar(name=args[i], **kw) for i in range(len(args))])
if lb is not None:
if np.isscalar(lb) or (isinstance(lb, np.ndarray) and lb.size == 1):
for v in r.view(np.ndarray):
v.lb = lb
else:
assert type(lb) in (list, tuple, ndarray)
for i, v in enumerate(r):
v.lb = lb[i]
if ub is not None:
if np.isscalar(ub) or (isinstance(ub, np.ndarray) and ub.size == 1):
for v in r.view(np.ndarray):
v.ub = ub
else:
assert type(ub) in (list, tuple, ndarray)
for i, v in enumerate(r):
v.ub = ub[i]
r._is_array_of_oovars = True
return r
def get_P(v, Lx, Ux, p, dataType, other=None, goal_is_nlh = True):
DefiniteRange = True
d = v.domain
if d is None:
raise FuncDesignerException('probably you are invoking boolean operation on continuous oovar')
if d is int or d is 'int':
raise FuncDesignerException('probably you are invoking boolean operation on non-boolean oovar')
inds = p._oovarsIndDict.get(v, None)
m = Lx.shape[0]
if inds is None:
# this oovar is fixed
res = {}
if v.domain is bool or v.domain is 'bool':
if goal_is_nlh:
T0 = True if p._x0[v] == 1 else False # 0 or 1
else:
T0 = False if p._x0[v] == 1 else True # 0 or 1
else:
assert other is not None, 'bug in FD kernel: called nlh with incorrect domain type'
if goal_is_nlh:
T0 = False if p._x0[v] != other else True
else:
T0 = False if p._x0[v] == other else True
return T0, res, DefiniteRange
#raise FuncDesignerException('probably you are trying to get nlh of fixed oovar, this is unimplemented in FD yet')
ind1, ind2 = inds
assert ind2-ind1 == 1, 'unimplemented for oovars of size > 1 yet'
lx, ux = Lx[:, ind1], Ux[:, ind1]
if d is bool or d is 'bool':
T0 = empty(m)
if goal_is_nlh:
T0.fill(inf)
T0[ux != lx] = 1.0 # lx = 0, ux = 1 => -log2(0.5) = 1
T0[lx == 1.0] = 0.0 # lx = 1 => ux = 1 => -log2(1) = 0
T2 = vstack((where(lx == 1, 0, inf), where(ux == 1, 0, inf))).T
else:
T0.fill(0)
T0[ux != lx] = 1.0
T0[lx == 1.0] = inf
T2 = vstack((where(lx == 1, inf, 0), where(ux == 1, inf, 0))).T
else:
assert other is not None, 'bug in FD kernel: called nlh with incorrect domain type'
mx = 0.5 * (lx + ux)
prev = 0
if prev:
ind = logical_and(mx==other, lx != ux)
if any(ind):
p.pWarn('seems like a categorical variables bug in FuncDesigner kernel, inform developers')
# mx[ind] += 1e-15 + 1e-15*abs(mx[ind])
I = searchsorted(d, lx, 'right') - 1
J = searchsorted(d, mx, 'right') - 1
#assert np.all(searchsorted(d, mx, 'right') == searchsorted(d, mx, 'left'))
K = searchsorted(d, ux, 'right') - 1
D0, D1, D2 = d[I], d[J], d[K]
d1, d2 = D0, D1
# if goal_is_nlh:
tmp1 = asfarray(J-I+1+where(d2==other, 1, 0))
tmp1[logical_or(otherd2)] = inf
# else:
# tmp1 = asfarray(J-I+where(d2==other, 0, 1))
# tmp1[logical_or(otherd2)] = 0
d1, d2 = D1, D2
# if goal_is_nlh:
tmp2 = asfarray(K-J+1+where(d2==other, 1, 0))
tmp2[logical_or(otherd2)] = inf
# else:
# tmp2 = asfarray(K-J+where(d2==other, 0, 1))
# tmp2[logical_or(otherd2)] = 0
if goal_is_nlh:
T2 = log2(vstack((tmp1, tmp2)).T)
else:
T2 = log2(vstack((tmp1, tmp2)).T)
d1, d2 = D0, D2
tmp = asfarray(K-I+where(d2==other, 1, 0))
tmp[logical_or(otherd2)] = inf
T0 = log2(tmp)
else:
assert np.all(d == array(d, int)) and len(d) == d[-1]-d[0]+1, 'bug in FD kernel'
#assert np.all(1e-6 < np.abs(np.array(mx, int)-mx))
assert goal_is_nlh, 'unimplemented yet'
tmp = ux - lx
tmp[other < lx] = inf
tmp[other > ux] = inf
tmp[logical_and(tmp==0, other == lx)] = 0.0
T0 = log2(tmp+1)
floor_mx = np.floor(mx)
tmp1 = floor_mx - lx
tmp1[other < lx] = inf
tmp1[other > floor_mx] = inf
tmp1[logical_and(tmp1==0, other == lx)] = 0.0
ceil_mx = np.ceil(mx)
tmp2 = ux - ceil_mx
tmp2[other > ux] = inf
tmp2[other < ceil_mx] = inf
tmp2[logical_and(tmp2==0, other == ux)] = 0.0
if goal_is_nlh:
T2 = log2(vstack((tmp1, tmp2)).T + 1.0)
else:
assert 0, 'unimplemented yet'
#T2 = log2(vstack((tmp1, tmp2)).T)
res = T2
return T0, res, DefiniteRange
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/overloads.py 000644 000766 000024 00000137545 12100431107 025774 0 ustar 00rfadams staff 000000 000000 PythonSum = sum
PythonMax = max
from ooFun import oofun
import numpy as np
from FDmisc import FuncDesignerException, Diag, Eye, raise_except, diagonal, DiagonalType
from ooFun import atleast_oofun, Vstack, Copy
from ooarray import ooarray
from Interval import TrigonometryCriticalPoints, nonnegative_interval, ZeroCriticalPointsInterval, box_1_interval
from numpy import atleast_1d, logical_and
from FuncDesigner.multiarray import multiarray
try:
from scipy.sparse import isspmatrix, lil_matrix as Zeros
scipyInstalled = True
except ImportError:
scipyInstalled = False
isspmatrix = lambda *args, **kw: False
Zeros = np.zeros
__all__ = []
#class unary_oofun_overload:
# def __init__(self, *args, **kwargs):
# assert len(args) == 1 and len(kwargs) == 0
# self.altFunc = args[0]
#
# def __call__(self, *args, **kwargs):
# assert len(args) == 1 and len(kwargs) == 0
# if not isinstance(args[0], oofun):
# return self.altFunc(*args)
# return self.fun(*args, **kwargs)
#@unary_oofun_overload(np.sin)
#def sin(inp):
# #if not isinstance(inp, oofun): return np.sin(inp)
# def d(x):
# return np.cos(x)
# return oofun(lambda x: np.sin(x), inp, d = d)
try:
import distribution
hasStochastic = True
except:
hasStochastic = False
#hasStochastic = False
st_sin = (lambda x: \
distribution.stochasticDistribution(sin(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([sin(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.sin(x))\
if hasStochastic\
else np.sin
def sin(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([sin(elem) for elem in inp])
elif hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(sin(inp.values), inp.probabilities.copy())._update(inp)
elif not isinstance(inp, oofun): return np.sin(inp)
return oofun(st_sin, inp,
d = lambda x: Diag(np.cos(x)),
vectorized = True,
criticalPoints = TrigonometryCriticalPoints)
#_interval = lambda domain: ufuncInterval(inp, domain, np.sin, TrigonometryCriticalPoints))
st_cos = (lambda x: \
distribution.stochasticDistribution(cos(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([cos(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.cos(x))\
if hasStochastic\
else np.cos
def cos(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([cos(elem) for elem in inp])
elif hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(cos(inp.values), inp.probabilities.copy())._update(inp)
elif not isinstance(inp, oofun): return np.cos(inp)
#return oofun(np.cos, inp, d = lambda x: Diag(-np.sin(x)))
return oofun(st_cos, inp,
d = lambda x: Diag(-np.sin(x)),
vectorized = True,
criticalPoints = TrigonometryCriticalPoints)
#_interval = lambda domain: ufuncInterval(inp, domain, np.cos, TrigonometryCriticalPoints))
st_tan = (lambda x: \
distribution.stochasticDistribution(tan(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([tan(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.tan(x))\
if hasStochastic\
else np.tan
def tan(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([tan(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(tan(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.tan(inp)
# TODO: move it outside of tan definition
def interval(*args):
raise 'interval for tan is unimplemented yet'
r = oofun(st_tan, inp, d = lambda x: Diag(1.0 / np.cos(x) ** 2), vectorized = True, interval = interval)
return r
__all__ += ['sin', 'cos', 'tan']
# TODO: cotan?
# TODO: rework it with matrix ops
get_box1_DefiniteRange = lambda lb, ub: logical_and(np.all(lb >= -1.0), np.all(ub <= 1.0))
st_arcsin = (lambda x: \
distribution.stochasticDistribution(arcsin(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([arcsin(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.arcsin(x))\
if hasStochastic\
else np.arcsin
def arcsin(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([arcsin(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(arcsin(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun):
return np.arcsin(inp)
r = oofun(st_arcsin, inp, d = lambda x: Diag(1.0 / np.sqrt(1.0 - x**2)), vectorized = True)
r.getDefiniteRange = get_box1_DefiniteRange
F_l, F_u = np.arcsin((-1, 1))
r._interval_ = lambda domain, dtype: box_1_interval(inp, np.arcsin, domain, dtype, F_l, F_u)
r.attach((inp>-1)('arcsin_domain_lower_bound_%d' % r._id, tol=-1e-7), (inp<1)('arcsin_domain_upper_bound_%d' % r._id, tol=-1e-7))
return r
st_arccos = (lambda x: \
distribution.stochasticDistribution(arccos(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([arccos(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.arccos(x))\
if hasStochastic\
else np.arccos
def arccos(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([arccos(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(arccos(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.arccos(inp)
r = oofun(st_arccos, inp, d = lambda x: Diag(-1.0 / np.sqrt(1.0 - x**2)), vectorized = True)
r.getDefiniteRange = get_box1_DefiniteRange
F_l, F_u = np.arccos((-1, 1))
r._interval_ = lambda domain, dtype: box_1_interval(inp, np.arccos, domain, dtype, F_l, F_u)
r.attach((inp>-1)('arccos_domain_lower_bound_%d' % r._id, tol=-1e-7), (inp<1)('arccos_domain_upper_bound_%d' % r._id, tol=-1e-7))
return r
st_arctan = (lambda x: \
distribution.stochasticDistribution(arctan(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([arctan(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.arctan(x))\
if hasStochastic\
else np.arctan
def arctan(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([arctan(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(arctan(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.arctan(inp)
return oofun(st_arctan, inp,
d = lambda x: Diag(1.0 / (1.0 + x**2)),
vectorized = True,
criticalPoints = False)
__all__ += ['arcsin', 'arccos', 'arctan']
st_sinh = (lambda x: \
distribution.stochasticDistribution(sinh(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([sinh(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.sinh(x))\
if hasStochastic\
else np.sinh
def sinh(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([sinh(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(sinh(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.sinh(inp)
return oofun(st_sinh, inp, d = lambda x: Diag(np.cosh(x)), vectorized = True, criticalPoints = False)
#def asdf(x):
## print (1, type(x))
## if isinstance(x, np.ndarray):
## print('1>', x.shape)
## print(x)
## if isinstance(x[0], np.ndarray):
## print(2, type(x[0]), x[0].shape)
## if isinstance(x, multiarray):
## print('-----')
## print (x.shape, x.view(np.ndarray).shape)
## print([type(elem) for elem in x.flat])
## print '===='
# r = distribution.stochasticDistribution(cosh(x.values), x.probabilities.copy())._update(x) \
# if isinstance(x, distribution.stochasticDistribution)\
# else np.array([cosh(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray)\
# else np.cosh(x)
# return r
#
#st_cosh = (asdf)\
#if hasStochastic\
#else np.cosh
st_cosh = \
(lambda x: \
distribution.stochasticDistribution(cosh(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([cosh(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.cosh(x))\
if hasStochastic\
else np.cosh
def cosh(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([cosh(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(cosh(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun):
return np.cosh(inp)
return oofun(st_cosh, inp, d = lambda x: Diag(np.sinh(x)), vectorized = True, _interval_=ZeroCriticalPointsInterval(inp, np.cosh))
__all__ += ['sinh', 'cosh']
st_tanh = (lambda x: \
distribution.stochasticDistribution(tanh(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([tanh(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.tanh(x))\
if hasStochastic\
else np.tanh
def tanh(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([tanh(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(tanh(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.tanh(inp)
return oofun(st_tanh, inp, d = lambda x: Diag(1.0/np.cosh(x)**2), vectorized = True, criticalPoints = False)
st_arctanh = (lambda x: \
distribution.stochasticDistribution(arctanh(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([arctanh(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.arctanh(x))\
if hasStochastic\
else np.arctanh
def arctanh(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([arctanh(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(arctanh(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.arctanh(inp)
r = oofun(st_arctanh, inp, d = lambda x: Diag(1.0/(1.0-x**2)), vectorized = True, criticalPoints = False)
r.getDefiniteRange = get_box1_DefiniteRange
r._interval_ = lambda domain, dtype: box_1_interval(inp, np.arctanh, domain, dtype, -np.inf, np.inf)
return r
__all__ += ['tanh', 'arctanh']
st_arcsinh = (lambda x: \
distribution.stochasticDistribution(arcsinh(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([arcsinh(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.arcsinh(x))\
if hasStochastic\
else np.arcsinh
def arcsinh(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([arcsinh(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(arcsinh(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.arcsinh(inp)
return oofun(st_arcsinh, inp, d = lambda x: Diag(1.0/np.sqrt(1+x**2)), vectorized = True, criticalPoints = False)
st_arccosh = (lambda x: \
distribution.stochasticDistribution(arccosh(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([arccosh(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.arccosh(x))\
if hasStochastic\
else np.arccosh
def arccosh(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([arccosh(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(arccosh(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.arccosh(inp)
r = oofun(st_arccosh, inp, d = lambda x: Diag(1.0/np.sqrt(x**2-1.0)), vectorized = True)
F0, shift = 0.0, 1.0
r._interval_ = lambda domain, dtype: nonnegative_interval(inp, np.arccosh, domain, dtype, F0, shift)
return r
__all__ += ['arcsinh', 'arccosh']
def angle(inp1, inp2):
# returns angle between 2 vectors
# TODO:
# 1) handle zero vector(s)
# 2) handle small numerical errors more properly
# (currently they are handled via constraint attached to arccos)
return arccos(sum(inp1*inp2)/sqrt(sum(inp1**2)*sum(inp2**2)))
st_exp = (lambda x: \
distribution.stochasticDistribution(exp(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([exp(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.exp(x))\
if hasStochastic\
else np.exp
def exp(inp):
if isinstance(inp, ooarray):
return ooarray([exp(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(exp(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.exp(inp)
return oofun(st_exp, inp, d = lambda x: Diag(np.exp(x)), vectorized = True, criticalPoints = False)
st_sqrt = (lambda x: \
distribution.stochasticDistribution(sqrt(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([sqrt(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.sqrt(x))\
if hasStochastic\
else np.sqrt
def sqrt(inp, attachConstraints = True):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([sqrt(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(sqrt(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun):
return np.sqrt(inp)
# def fff(x):
# print x
# return np.sqrt(x)
r = oofun(st_sqrt, inp, d = lambda x: Diag(0.5 / np.sqrt(x)), vectorized = True)
F0 = 0.0
r._interval_ = lambda domain, dtype: nonnegative_interval(inp, np.sqrt, domain, dtype, F0)
if attachConstraints: r.attach((inp>0)('sqrt_domain_zero_bound_%d' % r._id, tol=-1e-7))
return r
__all__ += ['angle', 'exp', 'sqrt']
st_abs = (lambda x: \
distribution.stochasticDistribution(abs(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([abs(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.abs(x))\
if hasStochastic\
else np.abs
def abs(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([abs(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(abs(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.abs(inp)
return oofun(st_abs, inp, d = lambda x: Diag(np.sign(x)), vectorized = True, _interval_ = ZeroCriticalPointsInterval(inp, np.abs))
#return oofun(np.abs, inp, d = lambda x: Diag(np.sign(x)), vectorized = True, criticalPoints = ZeroCriticalPoints)
__all__ += ['abs']
def log_interval(logfunc, inp):
def interval(domain, dtype):
lb_ub, definiteRange = inp._interval(domain, dtype)
lb, ub = lb_ub[0], lb_ub[1]
t_min, t_max = atleast_1d(logfunc(lb)), atleast_1d(logfunc(ub))
ind = lb < 0
if np.any(ind):
ind2 = ub>0
t_min[atleast_1d(logical_and(ind, ind2))] = -np.inf
definiteRange = False
ind = atleast_1d(ub == 0)
if np.any(ind):
t_max[ind] = np.nan
t_min[ind] = np.nan
definiteRange = False
#print definiteRange
# TODO: rework definiteRange with matrix operations
return np.vstack((t_min, t_max)), definiteRange
return interval
st_log = (lambda x: \
distribution.stochasticDistribution(log(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([log(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.log(x))\
if hasStochastic\
else np.log
def log(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([log(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(log(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.log(inp)
r = oofun(st_log, inp, d = lambda x: Diag(1.0/x), vectorized = True, _interval_ = log_interval(np.log, inp))
r.attach((inp>1e-300)('log_domain_zero_bound_%d' % r._id, tol=-1e-7))
return r
st_log10 = (lambda x: \
distribution.stochasticDistribution(log10(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([log10(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.log10(x))\
if hasStochastic\
else np.log10
INV_LOG_10 = 1.0 / np.log(10)
def log10(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([log10(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(log10(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.log10(inp)
r = oofun(st_log10, inp, d = lambda x: Diag(INV_LOG_10 / x), vectorized = True, _interval_ = log_interval(np.log10, inp))
r.attach((inp>1e-300)('log10_domain_zero_bound_%d' % r._id, tol=-1e-7))
return r
st_log2 = (lambda x: \
distribution.stochasticDistribution(log2(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([log2(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.log2(x))\
if hasStochastic\
else np.log2
INV_LOG_2 = 1.0 / np.log(2)
def log2(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([log2(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(log2(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun): return np.log2(inp)
r = oofun(st_log2, inp, d = lambda x: Diag(INV_LOG_2/x), vectorized = True, _interval_ = log_interval(np.log2, inp))
r.attach((inp>1e-300)('log2_domain_zero_bound_%d' % r._id, tol=-1e-7))
return r
__all__ += ['log', 'log2', 'log10']
def dot(inp1, inp2):
if not isinstance(inp1, oofun) and not isinstance(inp2, oofun): return np.dot(inp1, inp2)
def aux_d(x, y):
if y.size == 1:
r = np.empty_like(x)
r.fill(y)
return Diag(r)
else:
return y
r = oofun(lambda x, y: x * y if x.size == 1 or y.size == 1 else np.dot(x, y), [inp1, inp2], d=(lambda x, y: aux_d(x, y), lambda x, y: aux_d(y, x)))
r.getOrder = lambda *args, **kwargs: (inp1.getOrder(*args, **kwargs) if isinstance(inp1, oofun) else 0) + (inp2.getOrder(*args, **kwargs) if isinstance(inp2, oofun) else 0)
#r.isCostly = True
return r
def cross(a, b):
if not isinstance(a, oofun) and not isinstance(b, oofun): return np.cross(a, b)
def aux_d(x, y):
assert x.size == 3 and y.size == 3, 'currently FuncDesigner cross(x,y) is implemented for arrays of length 3 only'
return np.array([[0, -y[2], y[1]], [y[2], 0, -y[0]], [-y[1], y[0], 0]])
r = oofun(lambda x, y: np.cross(x, y), [a, b], d=(lambda x, y: -aux_d(x, y), lambda x, y: aux_d(y, x)))
r.getOrder = lambda *args, **kwargs: (a.getOrder(*args, **kwargs) if isinstance(a, oofun) else 0) + (b.getOrder(*args, **kwargs) if isinstance(b, oofun) else 0)
return r
__all__ += ['dot', 'cross']
def ceil(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([ceil(elem) for elem in inp])
if not isinstance(inp, oofun): return np.ceil(inp)
r = oofun(lambda x: np.ceil(x), inp, vectorized = True)
r._D = lambda *args, **kwargs: raise_except('derivative for FD ceil is unimplemented yet')
r.criticalPoints = False#lambda arg_infinum, arg_supremum: [np.ceil(arg_supremum)]
return r
def floor(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([floor(elem) for elem in inp])
if not isinstance(inp, oofun): return np.floor(inp)
r = oofun(lambda x: np.floor(x), inp, vectorized = True)
r._D = lambda *args, **kwargs: raise_except('derivative for FD floor is unimplemented yet')
r.criticalPoints = False#lambda arg_infinum, arg_supremum: [np.floor(arg_infinum)]
return r
st_sign = (lambda x: \
distribution.stochasticDistribution(sign(x.values), x.probabilities.copy())._update(x) \
if isinstance(x, distribution.stochasticDistribution)\
else np.array([sign(elem) for elem in x.flat]).view(multiarray) if isinstance(x, multiarray) and isinstance(x.flat[0], distribution.stochasticDistribution)
else np.sign(x))\
if hasStochastic\
else np.sign
def sign(inp):
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]):
return ooarray([sign(elem) for elem in inp])
if hasStochastic and isinstance(inp, distribution.stochasticDistribution):
return distribution.stochasticDistribution(sign(inp.values), inp.probabilities.copy())._update(inp)
if not isinstance(inp, oofun):
return np.sign(inp)
r = oofun(st_sign, inp, vectorized = True, d = lambda x: 0.0)
r.criticalPoints = False
return r
__all__ += ['ceil', 'floor', 'sign']
def sum_engine(r0, *args):
if not hasStochastic:
return PythonSum(args) + r0
Args, Args_st = [], {}
for elem in args:
if isinstance(elem, distribution.stochasticDistribution):
stDep = frozenset(elem.stochDep.keys())
tmp = Args_st.get(stDep, None)
if tmp is None:
Args_st[stDep] = [elem]
else:
Args_st[stDep].append(elem)
else:
Args.append(elem)
r = PythonSum(Args) + r0
if len(Args_st) == 0:
return r
# temporary
for key, val in Args_st.items():
maxDistributionSize = val[0].maxDistributionSize
break
stValues = Args_st.values()
# stValues = Args_st.values()
# T = list(set(stValues))[0]
# maxDistributionSize = next(iter(T)).maxDistributionSize
r1 = 0.0
for elem in stValues:
tmp = PythonSum(elem)
r1 = tmp + r1
r1.reduce(maxDistributionSize)
r1 = r1 + r
r1.maxDistributionSize = maxDistributionSize
return r1
def sum_interval(R0, r, INP, domain, dtype):
v = domain.modificationVar
if v is not None:
# self already must be in domain.storedSums
R, DefiniteRange = domain.storedSums[r][-1]
if not np.all(np.isfinite(R)):
R = np.asarray(R0, dtype).copy()
if domain.isMultiPoint:
R = np.tile(R, (1, len(list(domain.values())[0][0])))
#R = np.tile(R, (1, len(domain.values()[0][0])))
DefiniteRange = True
#####################
# !!! don't use sum([inp._interval(domain, dtype) for ...]) here
# to reduce memory consumption
for inp in INP:
arg_lb_ub, definiteRange = inp._interval(domain, dtype)
# DefiniteRange = logical_and(DefiniteRange, definiteRange)
if R.shape == arg_lb_ub.shape:
R += arg_lb_ub
else:
R = R + arg_lb_ub
#####################
return R, DefiniteRange
R=R.copy()
for inp in r.storedSumsFuncs[v]:
# TODO: mb rework definiteRange processing ?
arg_lb_ub, definiteRange = inp._interval(domain, dtype)
DefiniteRange = logical_and(DefiniteRange, definiteRange)
R += arg_lb_ub
R -= domain.storedSums[r][v]
# To supress inf-inf=nan, however, it doesn't work properly yet, other code is used
if np.any(np.isinf(arg_lb_ub)):
R[arg_lb_ub == np.inf] = np.inf
R[arg_lb_ub == -np.inf] = -np.inf
return R, DefiniteRange
else:
domain.storedSums[r] = {}
#assert np.asarray(r0).ndim <= 1
#R = np.asarray(R0, dtype).copy()
R = np.asarray(R0).copy()
if domain.isMultiPoint:
R = np.tile(R, (1, len(list(domain.values())[0][0])))
#####################
# !!! don't use sum([inp._interval(domain, dtype) for ...]) here
# to reduce memory consumption
DefiniteRange = True
# R_ = []
D = domain.storedSums[r]
for inp in INP:
arg_lb_ub, definiteRange = inp._interval(domain, dtype)
Tmp = inp._getDep() if not inp.is_oovar else [inp]
for oov in Tmp:
tmp = D.get(oov, None)
if tmp is None:
D[oov] = arg_lb_ub.copy()
else:
try:
D[oov] += arg_lb_ub
except:
# may be of different shape, e.g. for a fixed variable
D[oov] = D[oov] + arg_lb_ub
DefiniteRange = logical_and(DefiniteRange, definiteRange)
if R.shape == arg_lb_ub.shape:
R += arg_lb_ub
else:
R = R + arg_lb_ub
#####################
if v is None:
domain.storedSums[r][-1] = R, DefiniteRange
return R, DefiniteRange
def sum_derivative(r_, r0, INP, dep, point, fixedVarsScheduleID, Vars=None, fixedVars = None, useSparse = 'auto'):
# TODO: handle involvePrevData
# TODO: handle fixed vars
r = {}
isSP = hasattr(point, 'maxDistributionSize') and point.maxDistributionSize != 0
for elem in INP:
if not elem.is_oovar and (elem.input is None or len(elem.input)==0 or elem.input[0] is None):
continue # TODO: get rid if None, use [] instead
if elem.discrete: continue
# TODO: code cleanup
if elem.is_oovar:
if (fixedVars is not None and elem in fixedVars) or (Vars is not None and elem not in Vars): continue
sz = np.asarray(point[elem]).size
tmpres = Eye(sz) if not isinstance(point[elem], multiarray) else np.ones(sz).view(multiarray)
r_val = r.get(elem, None)
if isSP:
if r_val is not None:
r_val.append(tmpres)
else:
r[elem] = [tmpres]
else:
if r_val is not None:
if sz != 1 and isinstance(r_val, np.ndarray) and not isinstance(tmpres, np.ndarray): # i.e. tmpres is sparse matrix
tmpres = tmpres.toarray()
elif not np.isscalar(r_val) and not isinstance(r_val, np.ndarray) and isinstance(tmpres, np.ndarray):
r[elem] = r_val.toarray()
Tmp = tmpres.resolve(True) if isspmatrix(r[elem]) and type(tmpres) == DiagonalType else tmpres
try:
r[elem] += Tmp
except:
r[elem] = r[elem] + Tmp
else:
# TODO: check it for oovars with size > 1
r[elem] = tmpres
else:
tmp = elem._D(point, fixedVarsScheduleID, Vars, fixedVars, useSparse = useSparse)
for key, val in tmp.items():
r_val = r.get(key, None)
if isSP:
if r_val is not None:
r_val.append(val)
else:
r[key] = [val]
else:
if r_val is not None:
if not np.isscalar(val) and isinstance(r_val, np.ndarray) and not isinstance(val, np.ndarray): # i.e. tmpres is sparse matrix
val = val.toarray()
elif not np.isscalar(r_val) and not isinstance(r_val, np.ndarray) and isinstance(val, np.ndarray):
r[key] = r_val.toarray()
if isspmatrix(r_val) and type(val) == DiagonalType:
val = val.resolve(True)
elif isspmatrix(val) and type(r_val) == DiagonalType:
r[key] = r_val.resolve(True)
# TODO: rework it
try:
r[key] += val
except:
r[key] = r_val + val
else:
r[key] = Copy(val)
if isSP:
for key, val in r.items():
r[key] = sum_engine(0.0, *val)
if useSparse is False:
for key, val in r.items():
#if np.isscalar(val): val = np.asfarray(val)
if hasattr(val, 'toarray'):# and not isinstance(val, multiarray): # i.e. sparse matrix
r[key] = val.toarray()
if not isSP:
# TODO: rework it, don't recalculate each time
Size = np.asarray(r0).size
if Size == 1 and not point.isMultiPoint:
if r_._lastFuncVarsID == fixedVarsScheduleID:
if not np.isscalar(r_._f_val_prev):
Size = r_._f_val_prev.size
else:
Size = np.asarray(r_._getFuncCalcEngine(point, Vars = Vars, fixedVars = fixedVars, fixedVarsScheduleID = fixedVarsScheduleID)).size
if Size != 1 and not point.isMultiPoint:
for key, val in r.items():
if not isinstance(val, diagonal):
if np.isscalar(val) or np.prod(val.shape) <= 1:
tmp = np.empty((Size, 1))
tmp.fill(val if np.isscalar(val) else val.item())
r[key] = tmp
elif val.shape[0] != Size:
tmp = np.tile(val, (Size, 1))
r[key] = tmp
# elif np.asarray(val).size !=1:
# raise_except('incorrect size in FD sum kernel')
return r
def sum_getOrder(INP, *args, **kwargs):
orders = [0]+[inp.getOrder(*args, **kwargs) for inp in INP]
return np.max(orders)
def sum(inp, *args, **kwargs):
if type(inp) == np.ndarray and inp.dtype != object:
return np.sum(inp, *args, **kwargs)
if isinstance(inp, ooarray) and inp.dtype != object:
inp = inp.view(np.ndarray)
cond_ooarray = isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)])
if cond_ooarray and inp.size == 1:
return np.asscalar(inp).sum()
condIterableOfOOFuns = type(inp) in (list, tuple) or cond_ooarray
if not isinstance(inp, oofun) and not condIterableOfOOFuns:
return np.sum(inp, *args, **kwargs)
if isinstance(inp, ooarray) and any([isinstance(elem, oofun) for elem in atleast_1d(inp)]): inp = inp.tolist()
if condIterableOfOOFuns:
d, INP, r0 = [], [], 0.0
for elem in inp: # TODO: mb use reduce() or something like that
if not isinstance(elem, oofun):
# not '+=' because size can be changed from 1 to another value
r0 = r0 + np.asanyarray(elem) # so it doesn't work for different sizes
continue
INP.append(elem)
if len(INP) == 0:
return r0
r = oofun(lambda *args: sum_engine(r0, *args), INP, _isSum = True)
r._summation_elements = INP if np.isscalar(r0) and r0 == 0.0 else INP + [r0]
r.storedSumsFuncs = {}
for inp in INP:
Dep = [inp] if inp.is_oovar else inp._getDep()
for v in Dep:
if v not in r.storedSumsFuncs:
r.storedSumsFuncs[v] = set()
r.storedSumsFuncs[v].add(inp)
# TODO: check for fixed inputs
r.getOrder = lambda *args, **kw: sum_getOrder(INP, *args, **kw)
R0 = np.tile(r0, (2, 1))
r._interval_ = lambda *args, **kw: sum_interval(R0, r, INP, *args, **kw)
r.vectorized = True
r_dep = r._getDep()
r._D = lambda *args, **kw: sum_derivative(r, r0, INP, r_dep, *args, **kw)
# r.isCostly = True
return r
else:
return inp.sum(*args, **kwargs)#np.sum(inp, *args, **kwargs)
def prod(inp, *args, **kwargs):
if not isinstance(inp, oofun): return np.prod(inp, *args, **kwargs)
if len(args) != 0 or len(kwargs) != 0:
raise FuncDesignerException('oofun for prod(x, *args,**kwargs) is not implemented yet')
#r.getOrder = lambda *args, **kwargs: prod([(1 if not isinstance(inp, oofun) else inp.getOrder(*args, **kwargs)) for inp in self.input])
return inp.prod()
# Todo: implement norm_1, norm_inf etc
def norm(*args, **kwargs):
if len(kwargs) or len(args) > 1:
return np.linalg.norm(*args, **kwargs)
r = sqrt(sum(args[0]**2), attachConstraints=False)
if isinstance(r, oofun):
r.hasDefiniteRange=True
return r
__all__ += ['sum', 'prod', 'norm']
#def stack(*args, **kwargs):
# assert len(kwargs) == 0 and len(args) != 0
# if len(args) == 1:
# assert type(args[0]) in (list, tuple)
# if not any([isinstance(arg, oofun) for arg in args[0]]): return np.hstack(args)
# #f = lambda *Args: np.hstack([arg(Args) if isinstance(arg, oofun) else arg for arg in args[0]])
# def f(*Args):
# r = np.hstack([arg.fun(Args) if isinstance(arg, oofun) else arg for arg in args[0]])
# print '1:', r
# raise 0
# return r
# #d = lambda *Args: np.hstack([arg.d(Args).reshape(-1, 1) if isinstance(arg, oofun) else np.zeros((len(arg))) for arg in args[0]])
# def d(*Args):
# r = np.hstack([arg.d(Args).reshape(-1, 1) if isinstance(arg, oofun) else np.zeros((len(arg))) for arg in args[0]])
# print '2:', r
# return r
# print 'asdf', args[0]
# return oofun(f, args[0], d=d)
# else:
# raise FuncDesignerException('unimplemented yet')
# #assert isinstance(args[0], oofun)
#def norm(inp, *args, **kwargs):
# if len(args) != 0 or len(kwargs) != 0:
# raise FuncDesignerException('oofun for norm(x, *args,**kwargs) is not implemented yet')
#
# #np.linalg.norm
# f = lambda x: np.sqrt(np.sum(x**2))
#
# r = oofun(f, inp, isCostly = True)
#
# def d(x):
#
#
# #r.d = lambda *args, **kwargs:
#
# #s = r(x)
# #return Diag(x / s if s != 0 else np.zeros(x.size)) # however, dirivative doesn't exist in (0,0,..., 0)
# r.d = d
#
# return r
def size(inp, *args, **kwargs):
if not isinstance(inp, oofun): return np.size(inp, *args, **kwargs)
return inp.size
def ifThenElse(condition, val1, val2, *args, **kwargs):
# for future implementation
assert len(args) == 0 and len(kwargs) == 0
Val1 = atleast_oofun(val1)#fixed_oofun(val1) if not isinstance(val1, oofun) else val1
#if np.isscalar(val1): raise 0
Val2 = atleast_oofun(val2)#fixed_oofun(val2) if not isinstance(val2, oofun) else val2
if isinstance(condition, bool):
return Val1 if condition else Val2
elif isinstance(condition, oofun):
f = lambda conditionResult, value1Result, value2Result: value1Result if conditionResult else value2Result
# !!! Don't modify it elseware function will evaluate both expressions despite of condition value
r = oofun(f, [condition, val1, val2])
r.D = lambda point, *args, **kwargs: (Val1.D(point, *args, **kwargs) if isinstance(Val1, oofun) else {}) if condition(point) else \
(Val2.D(point, *args, **kwargs) if isinstance(Val2, oofun) else {})
r._D = lambda point, *args, **kwargs: (Val1._D(point, *args, **kwargs) if isinstance(Val1, oofun) else {}) if condition(point) else \
(Val2._D(point, *args, **kwargs) if isinstance(Val2, oofun) else {})
r.d = errFunc
# TODO: try to set correct value from val1, val2 if condition is fixed
# def getOrder(Vars=None, fixedVars=None, *args, **kwargs):
# dep = condition.getDep()
# if Vars is not None and dep.is
return r
else:
raise FuncDesignerException('ifThenElse requires 1st argument (condition) to be either boolean or oofun, got %s instead' % type(condition))
__all__ += ['size', 'ifThenElse']
def decision(*args, **kwargs):
pass
def max(inp, *args, **kwargs):
if type(inp) in (list, tuple, np.ndarray) \
and (len(args) == 0 or len(args) == 1 and not isinstance(args[0], oofun)) \
and not any([isinstance(elem, oofun) for elem in (inp if type(inp) in (list, tuple) else np.atleast_1d(inp))]):
return np.max(inp, *args, **kwargs)
assert len(args) == len(kwargs) == 0, 'incorrect data type in FuncDesigner max or not implemented yet'
if isinstance(inp, oofun):
f = lambda x: np.max(x)
# def f(x):
# print np.max(x)
# return np.max(x)
def d(x):
df = inp.d(x)
ind = np.argmax(x)
return df[ind, :]
def interval(domain, dtype):
lb_ub, definiteRange = inp._interval(domain, dtype)
tmp1, tmp2 = lb_ub[0], lb_ub[1]
return np.vstack((np.max(np.vstack(tmp1), 0), np.max(np.vstack(tmp2), 0))), np.all(definiteRange, 0)
r = oofun(f, inp, d = d, size = 1, _interval_ = interval)
elif type(inp) in (list, tuple, ooarray):
f = lambda *args: np.max([arg for arg in args])
def interval(domain, dtype):
arg_inf, arg_sup, tmp, DefiniteRange = [], [], -np.inf, True
for _inp in inp:
if isinstance(_inp, oofun):
#tmp1, tmp2 = _inp._interval(domain, dtype)
lb_ub, definiteRange = _inp._interval(domain, dtype)
tmp1, tmp2 = lb_ub[0], lb_ub[1]
arg_inf.append(tmp1)
arg_sup.append(tmp2)
DefiniteRange = logical_and(DefiniteRange, definiteRange)
elif tmp < _inp:
tmp = _inp
r1, r2 = np.max(np.vstack(arg_inf), 0), np.max(np.vstack(arg_sup), 0)
r1[r1 _inp:
tmp = _inp
r1, r2 = np.min(np.vstack(arg_inf), 0), np.min(np.vstack(arg_sup), 0)
if np.isfinite(tmp):
r1[r1>tmp] = tmp
r2[r2>tmp] = tmp
return np.vstack((r1, r2)), DefiniteRange
r = oofun(f, inp, size = 1, _interval_ = interval)
def _D(point, *args, **kwargs):
ind = np.argmin([(s(point) if isinstance(s, oofun) else s) for s in r.input])
return r.input[ind]._D(point, *args, **kwargs) if isinstance(r.input[ind], oofun) else {}
r._D = _D
else:
return np.min(inp, *args, **kwargs)
return r
__all__ += ['min', 'max']
#def fixed_oofun(Val):
# val = np.asfarray(Val)
# f = lambda: Val
# r = oofun(f, input=[])
# r._D = lambda *args, **kwargs: {}
# r.D = lambda *args, **kwargs: {}
# r.discrete = True
# return r
det3 = lambda a, b, c: a[0] * (b[1]*c[2] - b[2]*c[1]) - a[1] * (b[0]*c[2] - b[2]*c[0]) + a[2] * (b[0]*c[1] - b[1]*c[0])
__all__ += ['det3']
def hstack(tup): # overload for oofun[ind]
c = [isinstance(t, (oofun, ooarray)) for t in tup]
if any([isinstance(t, ooarray) for t in tup]):
return ooarray(np.hstack(tup))
if not any(c):
return np.hstack(tup)
#an_oofun_ind = np.where(c)[0][0]
f = lambda *x: np.hstack(x)
# def d(*x):
#
# r = [elem.d(x[i]) if c[i] else None for i, elem in enumerate(tup)]
# size = atleast_1d(r[an_oofun_ind]).shape[0]
# r2 = [elem if c[i] else Zeros(size) for elem in r]
# return r2
#= lambda *x: np.hstack([elem.d(x) if c[i] else elem for elem in tup])
# f = lambda x: x[ind]
# def d(x):
# Xsize = Len(x)
# condBigMatrix = Xsize > 100
# if condBigMatrix and scipyInstalled:
# r = SparseMatrixConstructor((1, x.shape[0]))
# r[0, ind] = 1.0
# else:
# if condBigMatrix and not scipyInstalled: self.pWarn(scipyAbsentMsg)
# r = zeros_like(x)
# r[ind] = 1
# return r
def getOrder(*args, **kwargs):
orders = [0]+[inp.getOrder(*args, **kwargs) for inp in tup]
return np.max(orders)
r = oofun(f, tup, getOrder = getOrder)
#!!!!!!!!!!!!!!!!! TODO: sparse
def _D(*args, **kwargs):
# TODO: rework it, especially if sizes are fixed and known
# TODO: get rid of fixedVarsScheduleID
sizes = [(t(args[0], fixedVarsScheduleID = kwargs.get('fixedVarsScheduleID', -1)) if c[i] else np.asarray(t)).size for i, t in enumerate(tup)]
tmp = [elem._D(*args, **kwargs) if c[i] else None for i, elem in enumerate(tup)]
res = {}
for v in r._getDep():
Temp = []
for i, t in enumerate(tup):
if c[i]:
temp = tmp[i].get(v, None)
if temp is not None:
Temp.append(temp if type(temp) != DiagonalType else temp.resolve(kwargs['useSparse']))
else:
# T = next(iter(tmp[i].values()))
# sz = T.shape[0] if type(T) == DiagonalType else np.atleast_1d(T).shape[0]
Temp.append((Zeros if sizes[i] * np.asarray(args[0][v]).size > 1000 else np.zeros)((sizes[i], np.asarray(args[0][v]).size)))
else:
sz = np.atleast_1d(t).shape[0]
Temp.append(Zeros((sz, 1)) if sz > 100 else np.zeros(sz))
rr = Vstack([elem for elem in Temp])
#print type(rr)
res[v] = rr if not isspmatrix(rr) or 0.3 * prod(rr.shape) > rr.size else rr.toarray()
#print type(res[v])
return res
r._D = _D
return r
__all__ += ['hstack']
# TODO: move the func into fdmisc.py
def errFunc(*args, **kwargs):
# this function shouldn't be ever called, an FD kernel hack has been involved
raise FuncDesignerException('error in FuncDesigner kernel, inform developers')
#for func in (sin, arctan):
# i0 = func._interval
# def f2(domain, dtype):
# if type(domain) == dict:
# return i0(domain, dtype)
# r = domain.storedIntervals.get(self, None)
# if r is None:
# r = i0(domain, dtype)
# domain.storedIntervals[self] = r
# return r
# func._interval = f2
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/sle.py 000644 000766 000024 00000006725 12100431107 024554 0 ustar 00rfadams staff 000000 000000 from FDmisc import FuncDesignerException
from ooFun import oofun, BaseFDConstraint
from numpy import nan, zeros, isscalar, inf
class sle:
# System of linear equations
_isInitialized = False
matrixSLEsolver = 'autoselect'
def __init__(self, equations, *args, **kwargsForOpenOptSLEconstructor):
if len(args) > 0: FuncDesignerException('incorrect sle definition, too many args are obtained')
if type(equations) not in [list, tuple, set]:
raise FuncDesignerException('argument of sle constructor should be Python tuple or list of equations or oofuns')
self.equations = equations
try:
from openopt import SLE
except:
s = "Currently to solve SLEs via FuncDesigner you should have OpenOpt installed; maybe in future the dependence will be ceased"
raise FuncDesignerException(s)
self.decodeArgs(*args, **kwargsForOpenOptSLEconstructor)
if 'iprint' not in kwargsForOpenOptSLEconstructor.keys():
kwargsForOpenOptSLEconstructor['iprint'] = -1
self.p = SLE(self.equations, self.startPoint, **kwargsForOpenOptSLEconstructor)
self.p._Prepare()
self.A, self.b = self.p.C, self.p.d
self.n = self.p.C.shape[0]
self.decode = lambda x: self.p._vector2point(x)
def solve(self, *args): # mb for future implementation - add **kwargsForOpenOptSLEconstructor here as well
if len(args) > 2:
raise FuncDesignerException('incorrect number of args, should be at most 2 (startPoint and/or solver name, order: any)')
self.decodeArgs(*args)
r = self.p.solve(matrixSLEsolver=self.matrixSLEsolver)
if r.istop >= 0:
return r
else:
R = {}
for key, value in self.p.x0.items():
R[key] = value * nan
r.xf = R
r.ff = inf
return r
def decodeArgs(self, *args, **kwargs):
hasStartPoint = False
for arg in args:
if isinstance(arg, str):
self.matrixSLEsolver = arg
elif isinstance(arg, dict):
startPoint = args[0]
hasStartPoint = True
else:
raise FuncDesignerException('incorrect arg type, should be string (solver name) or dict (start point)')
if 'startPoint' in kwargs:
startPoint = kwargs['startPoint']
hasStartPoint = True
if not hasStartPoint:
if hasattr(self, 'startPoint'): return # established from __init__
involvedOOVars = set()
for Elem in self.equations:
elem = Elem.oofun if Elem.isConstraint else Elem
if elem.is_oovar:
involvedOOVars.add(elem)
else:
involvedOOVars.update(elem._getDep())
startPoint = {}
for oov in involvedOOVars:
if isscalar(oov.size):
startPoint[oov] = zeros(oov.size)
else:
startPoint[oov] = 0
self.startPoint = startPoint
def eig(self, *args, **kw):
from openopt import EIG
p = EIG(self.A, *args, **kw)
r = p.solve()
vectors = []
for i in range(len(p.eigenvalues)):
vectors.append(self.decode(p.eigenvectors[:, i]))
r.eigenvectors = vectors
return r
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/stencils.py 000644 000766 000024 00000004051 12100431107 025603 0 ustar 00rfadams staff 000000 000000 from numpy import all, abs
from overloads import hstack
from FDmisc import FuncDesignerException
def d(arg, v, **kw):#, *args, **kw):
N = len(v)
# if len(args) == 1:
# derivativeSide = args[0]
# assert derivativeSide in ('left', 'right', 'both')
# else:
# derivativeSide = 'both'
stencil = kw.get('stencil', 3)
if stencil not in (2, 3):
raise FuncDesignerException('for d1 only stencil = 2 and 3 are implemented')
timestep = v[1]-v[0]
if not all(abs(v[1:] - v [:-1] - timestep) < 1e-10):
raise FuncDesignerException('unimplemented for non-uniform step yet')
if stencil == 2:
r1 = -3*arg[0] + 4*arg[1] - arg[2]
r2 = (arg[2:N] - arg[0:N-2]) / 2.0
r3 = 3*arg[N-1] - 4*arg[N-2] + arg[N-3]
return hstack((r1, r2, r3)) / timestep
elif stencil == 3:
r1 = -22 * arg[0] + 36 * arg[1] - 18 * arg[2] + 4 * arg[3]
r2 = -22 * arg[1] + 36 * arg[2] - 18 * arg[3] + 4 * arg[4] # TODO: mb rework it?
r3 = arg[0:N-4] -8*arg[1:N-3] + 8*arg[3:N-1] - arg[4:N]
r4 = 22 * arg[N-5] - 36 * arg[N-4] + 18 * arg[N-3] - 4 * arg[N-2] # TODO: mb rework it?
r5 = 22 * arg[N-4] - 36 * arg[N-3] + 18 * arg[N-2] - 4 * arg[N-1]
return hstack((r1, r2, r3, r4, r5)) / (12*timestep)
# if derivativeSide == 'both':
# r = hstack((r1, r2, r3))
# elif derivativeSide == 'left':
# r = hstack((r1, r2))
# else: # derivativeSide == 'right'
# r = hstack((r2, r3))
return r
def d2(arg, v, **kw):#, *args, **kw):
N = len(v)
timestep = v[1]-v[0]
if not all(abs(v[1:] - v [:-1] - timestep) < 1e-10):
raise FuncDesignerException('unimplemented for non-uniform step yet')
stencil = kw.get('stencil', 1)
if stencil not in (1, ):
raise FuncDesignerException('for d2 only stencil = 1 is implemented')
if stencil == 1:
r1 = arg[0] - 2*arg[1] + arg[2]
r2 = arg[0:N-2] - 2 * arg[1:N-1] + arg[2:N]
r3 = arg[N-1] - 2*arg[N-2] + arg[N-3]
return hstack((r1, r2, r3)) / timestep**2
vendor/oosuite/FuncDesigner/build/lib/FuncDesigner/translator.py 000644 000766 000024 00000017763 12100431107 026166 0 ustar 00rfadams staff 000000 000000 # Handling of FuncDesigner probs
from numpy import empty, hstack, vstack, asfarray, all, atleast_1d, cumsum, asarray, zeros, ndarray, prod, ones, isscalar, nan, array_equal, copy, array
#from nonOptMisc import scipyInstalled, Hstack, Vstack, Find, isspmatrix, SparseMatrixConstructor, DenseMatrixConstructor, Bmat
from FDmisc import FuncDesignerException
from ooPoint import ooPoint
DenseMatrixConstructor = lambda shape: zeros(shape)
isspmatrix = lambda *args: False
class FuncDesignerTranslator:
# freeVars = []
# fixedVars = []
def __init__(self, PointOrVariables, **kwargs): #, freeVars=None, fixedVars=None
#assert freeVars is not None or fixedVars is not None, 'at most one parameter of "fixedVars" and "freeVars" is allowed'
#assert 'freeVars' not in kwargs, 'only "fixedVars" and "freeVars" arguments are allowed, not "freeVars"'
if isinstance(PointOrVariables, dict):
Point = PointOrVariables
Variables = Point.keys()
self._sizeDict = dict([(v, asarray(PointOrVariables[v]).size) for v in PointOrVariables])
self._shapeDict = dict([(v, asarray(PointOrVariables[v]).shape) for v in PointOrVariables])
# TODO: assert v.size (if provided) == PointOrVariables[v]).size
# and same with shapes
else:
assert type(PointOrVariables) in [list, tuple, set]
Variables = PointOrVariables
self._sizeDict = dict([(v, (v.size if hasattr(v, 'size') and isinstance(v.size, int) else 1)) for v in Variables])
self._shapeDict = dict([(v, (v.shape if hasattr(v, 'shape') else ())) for v in Variables])
self._variables = Variables
self.n = sum(self._sizeDict.values())
oovar_sizes = list(self._sizeDict.values()) # FD: for opt oovars only
oovar_indexes = cumsum([0] + oovar_sizes)
self.oovarsIndDict = dict([(v, (oovar_indexes[i], oovar_indexes[i+1])) for i, v in enumerate(Variables)])
# TODO: mb use oovarsIndDict here as well (as for derivatives?)
#startDictData = [] #if fixedVars is None else [(v, startPoint[v]) for v in fixedVars]
# TODO: involve fixed variables
self._SavedValues = {'prevX':nan}
def vector2point(x):
isComplexArray = isinstance(x, ndarray) and str(x.dtype).startswith('complex')
if isComplexArray:
x = atleast_1d(array(x, copy=True))
else:
x = atleast_1d(array(x, copy=True, dtype=float))
if array_equal(x, self._SavedValues['prevX']):
return self._SavedValues['prevVal']
# without copy() ipopt and probably others can replace it by noise after closing
kw = {'skipArrayCast':True} if isComplexArray else {}
r = ooPoint([(v, x[oovar_indexes[i]:oovar_indexes[i+1]]) for i, v in enumerate(self._variables)], **kw)
self._SavedValues['prevVal'] = r
self._SavedValues['prevX'] = copy(x)
return r
self.vector2point = vector2point
point2vector = lambda self, point: array(atleast_1d(hstack([(point[v] if v in point else zeros(self._shapeDict[v])) for v in self._variables])), 'float')
def pointDerivative2array(self, pointDerivarive, useSparse = False, func=None, point=None):
# useSparse can be True, False, 'auto'
# !!!!!!!!!!! TODO: implement useSparse = 'auto' properly
if useSparse == 'auto' and not scipyInstalled:
useSparse = False
if useSparse is not False and not scipyInstalled:
raise FuncDesignerException('to handle sparse matrices you should have module scipy installed')
# however, this check is performed in other function (before this one)
# and those constraints are excluded automaticvally
n = self.n
if len(pointDerivarive) == 0:
if func is not None:
assert point is not None
funcLen = func(point).size
if useSparse:
return SparseMatrixConstructor((funcLen, n))
else:
return DenseMatrixConstructor((funcLen, n))
else:
raise FuncDesignerException('unclear error, maybe you have constraint independend on any optimization variables')
key, val = list(pointDerivarive.items())[0]
if isscalar(val) or (isinstance(val, ndarray) and val.shape == ()):
val = atleast_1d(val)
var_inds = self.oovarsIndDict[key]
# val.size works in other way (as nnz) for scipy.sparse matrices
funcLen = int(round(prod(val.shape) / (var_inds[1] - var_inds[0])))
newStyle = 1
# TODO: remove "useSparse = False", replace by code from FDmisc
assert useSparse is False, 'sparsity is not implemented in FD translator yet'
if useSparse is not False and newStyle:
r2 = []
hasSparse = False
zeros_start_ind = 0
zeros_end_ind = 0
for i, var in enumerate(freeVars):
if var in pointDerivarive:#i.e. one of its keys
if zeros_end_ind != zeros_start_ind:
r2.append(SparseMatrixConstructor((funcLen, zeros_end_ind - zeros_start_ind)))
zeros_start_ind = zeros_end_ind
tmp = pointDerivarive[var]
if isspmatrix(tmp):
hasSparse = True
else:
tmp = asarray(tmp) # else bug with scipy sparse hstack
if tmp.ndim < 2:
tmp = tmp.reshape(funcLen, prod(tmp.shape) // funcLen)
r2.append(tmp)
else:
zeros_end_ind += oovar_sizes[i]
hasSparse = True
if zeros_end_ind != zeros_start_ind:
r2.append(SparseMatrixConstructor((funcLen, zeros_end_ind - zeros_start_ind)))
r3 = Hstack(r2) if hasSparse else hstack(r2)
if isspmatrix(r3) and r3.nnz > 0.25 * prod(r3.shape): r3 = r3.A
return r3
# r2 = []
# hasSparse = False
# for i, var in enumerate(freeVars):
# if var in pointDerivarive:#i.e. one of its keys
# tmp = pointDerivarive[var]
# if isspmatrix(tmp): hasSparse = True
# if isinstance(tmp, float) or (isinstance(tmp, ndarray) and tmp.shape == ()):
# tmp = atleast_1d(tmp)
# if tmp.ndim < 2:
# tmp = tmp.reshape(funcLen, prod(tmp.shape) // funcLen)
# r2.append(tmp)
# else:
# r2.append(SparseMatrixConstructor((funcLen, oovar_sizes[i])))
# hasSparse = True
# r3 = Hstack(r2) if hasSparse else hstack(r2)
# if isspmatrix(r3) and r3.nnz > 0.25 * prod(r3.shape): r3 = r3.A
# return r3
else:
if funcLen == 1:
r = DenseMatrixConstructor(n)
else:
if useSparse:
r = SparseMatrixConstructor((n, funcLen))
else:
r = DenseMatrixConstructor((n, funcLen))
for key, val in pointDerivarive.items():
# TODO: remove indexes, do as above for sparse
indexes = self.oovarsIndDict[key]
if not useSparse and isspmatrix(val): val = val.A
if r.ndim == 1:
r[indexes[0]:indexes[1]] = val if isscalar(val) else val.flatten()
else:
r[indexes[0]:indexes[1], :] = val.T
if useSparse is True and funcLen == 1:
return SparseMatrixConstructor(r)
else:
return r.T if r.ndim > 1 else r.reshape(1, -1)
vendor/oosuite/._INSTALL.txt 000755 000766 000024 00000000324 12100431107 016676 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/INSTALL.txt 000755 000766 000024 00000001663 12100431107 016470 0 ustar 00rfadams staff 000000 000000 OpenOpt, DerApproximator, FuncDesigner are packages written in Python language (www.python.org) + numpy (www.scipy.org/numpy).
They are stand-alone modules, but DerApproximator is required by some nonlinear solvers from OpenOpt and during automatic differentiation for those FuncDesigner oofuns w/o supplied routines for obtaining derivatives, and FuncDesigner models can be solved/optimized by OpenOpt.
To install them all you could just run (with admin rights, + you should have package "python-setuptools" installed)
python install_all
or
python develop_all
Former goes to directories OpenOpt, FuncDesigner, DerApproximator and invokes "python setup.py install" in each of them, latter do same the things with "python setup.py develop" (in-place installation).
For involving installation options beyond defaults you should read setuptools doc, go to these directories and perform "python setup.py ..." with desired options.
HTH,
Dmitrey.
vendor/oosuite/._OpenOpt 000755 000766 000024 00000000324 12100432771 016346 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/ 000755 000766 000024 00000000000 12100432771 016204 5 ustar 00rfadams staff 000000 000000 vendor/oosuite/OpenOpt/build/ 000755 000766 000024 00000000000 12100432771 017303 5 ustar 00rfadams staff 000000 000000 vendor/oosuite/OpenOpt/._COPYING.txt 000755 000766 000024 00000000324 12100431107 020264 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/COPYING.txt 000755 000766 000024 00000001773 12100431107 020060 0 ustar 00rfadams staff 000000 000000 Copyright (c) 2007-2009, Dmitrey Kroshko, www.icyb.kiev.ua optimization department
OpenOpt Kernel license is New BSD license, allows using from both open- and closed-code software
All connected solvers have their own licenses.
-----------------------------------
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
vendor/oosuite/OpenOpt/dist/ 000755 000766 000024 00000000000 12100432771 017147 5 ustar 00rfadams staff 000000 000000 vendor/oosuite/OpenOpt/._DOCUMENTATION.html 000755 000766 000024 00000000324 12100431107 021312 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/DOCUMENTATION.html 000755 000766 000024 00000000110 12100431107 021066 0 ustar 00rfadams staff 000000 000000
vendor/oosuite/OpenOpt/._INSTALL.txt 000755 000766 000024 00000000324 12100431107 020262 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/INSTALL.txt 000755 000766 000024 00000001676 12100431107 020060 0 ustar 00rfadams staff 000000 000000 Do you have any problems with installation for Linux, MS Windows, MacOS?
Check webpage
http://openopt.org/Install
Installation (very briefly, link above leads to more precize):
1. Ensure you have Python (v >= 2.5 is recommended), numpy (v >= 1.1.0 is recommended) and setuptools installed.
You can use Linux software channels, like this for Debian/*UBUNTU:
[sudo] aptitude install python-numpy
[sudo] aptitude install python-setuptools
(Linux OSes already have Python installed)
# optional, for graphics output:
[sudo] aptitude install python-matplotlib
# optional packages that may be present in software channels:
# python-scipy, python-cvxopt, some others
All optional packages and/or solvers could be installed after OO installation, when they will be required for the tasks involved.
2.
[sudo] python setup.py install
or
[sudo] python setup.py develop
3. To check installation you could try running "python nlp_1.py" from /examples directory
vendor/oosuite/OpenOpt/._openopt 000755 000766 000024 00000000324 12100432125 020023 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/openopt/ 000755 000766 000024 00000000000 12100432125 017661 5 ustar 00rfadams staff 000000 000000 vendor/oosuite/OpenOpt/openopt.egg-info/ 000755 000766 000024 00000000000 12100432771 021362 5 ustar 00rfadams staff 000000 000000 vendor/oosuite/OpenOpt/._README.txt 000755 000766 000024 00000000324 12100431107 020111 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/README.txt 000755 000766 000024 00000000721 12100431107 017675 0 ustar 00rfadams staff 000000 000000 openopt license: BSD
openopt usage:
from openopt import *
#and then see
help(NLP), help(LP), help(QP) # ... etc
See also directory "examples", using on-line documentation is highly recommended:
http://openopt.org/Doc
See also:
OpenOpt homepage:
http://openopt.org
New! Numerical optimization forum:
http://forum.openopt.org/
Found OpenOpt useful?
Please mention it in our guestbook
http://forum.openopt.org/viewforum.php?id=11
Regards,
OpenOpt developers
vendor/oosuite/OpenOpt/._setup.cfg 000755 000766 000024 00000000324 12100431107 020234 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/setup.cfg 000755 000766 000024 00000000064 12100431107 020020 0 ustar 00rfadams staff 000000 000000 [egg_info]
#tag_svn_revision = 1
#tag_build = .dev
vendor/oosuite/OpenOpt/._setup.py 000755 000766 000024 00000000324 12100431107 020125 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/setup.py 000755 000766 000024 00000007375 12100431107 017725 0 ustar 00rfadams staff 000000 000000 #! /usr/bin/env python
descr = """
"""
from os.path import join
import os
import sys, compileall
DISTNAME = 'openopt'
DESCRIPTION = 'A python module for numerical optimization'
LONG_DESCRIPTION = descr
MAINTAINER = 'Dmitrey Kroshko',
MAINTAINER_EMAIL = 'dmitrey-at-openopt-dot-org',
URL = 'http://openopt.org',
LICENSE = 'new BSD'
sys.path.append(os.getcwd() + os.sep + 'openopt')
from ooVersionNumber import __version__ as ooVer
openopt_version = ooVer
DOWNLOAD_URL = 'http://openopt.org/images/3/33/OpenOpt.zip'
try:
import setuptools
except:
print('you should have setuptools installed (http://pypi.python.org/pypi/setuptools), for some Linux distribs you can get it via [sudo] apt-get install python-setuptools')
print('press Enter for exit...')
raw_input()
exit()
from distutils.errors import DistutilsError
#from numpy.distutils.system_info import system_info, NotFoundError, dict_append, so_ext
from numpy.distutils.core import setup, Extension
DOC_FILES = []
def configuration(parent_package='',top_path=None, package_name=DISTNAME):
if os.path.exists('MANIFEST'): os.remove('MANIFEST')
pkg_prefix_dir = 'openopt'
# Get the version
from numpy.distutils.misc_util import Configuration
config = Configuration(package_name,parent_package,top_path,
version = openopt_version,
maintainer = MAINTAINER,
maintainer_email = MAINTAINER_EMAIL,
description = DESCRIPTION,
license = LICENSE,
url = URL,
download_url = DOWNLOAD_URL,
long_description = LONG_DESCRIPTION)
# XXX: once in SVN, should add svn version...
#print config.make_svn_version_py()
# package_data does not work with sdist for setuptools 0.5 (setuptools bug),
# so we need to add them here while the bug is not solved...
return config
if __name__ == "__main__":
# solverPaths = {}
# #File = string.join(__file__.split(os.sep)[:-1], os.sep)
# for root, dirs, files in os.walk('openopt'+os.sep +'solvers'):
# #for root, dirs, files in os.walk(os.path.dirname(file)+os.sep+'solvers'):
# rd = root.split(os.sep)
# if '.svn' in rd: continue
# rd = rd[rd.index('solvers')+1:]
# for file in files:
# if len(file)>6 and file[-6:] == '_oo.py':
# solverPaths[file[:-6]] = string.join(rd,'.') + '.'+file[:-3]
# f = open('solverPaths.py', 'w')
# f.write('solverPaths = ' + str(solverPaths))
# f.close()
# shutil.move('solverPaths.py', 'openopt' + os.sep + 'kernel' + os.sep + 'solverPaths.py')
# setuptools version of config script
# package_data does not work with sdist for setuptools 0.5 (setuptools bug)
# So we cannot add data files via setuptools yet.
#data_files = ['test_data/' + i for i in TEST_DATA_FILES]
#data_files.extend(['docs/' + i for i in doc_files])
setup(configuration = configuration,
install_requires=['numpy','setproctitle'], # can also add version specifiers #namespace_packages=['kernel'],
#setup_requires = 'setproctitle',
#py_modules = ['kernel', 'tests', 'examples', 'solvers'],
packages=setuptools.find_packages(),
include_package_data = True,
#package_data = '*.txt',
test_suite='',#"openopt.tests", # for python setup.py test
zip_safe=False, # the package can run out of an .egg file
#FIXME url, download_url, ext_modules
classifiers =
[ 'Development Status :: 4 - Beta',
'Environment :: Console',
'Intended Audience :: Developers',
'Intended Audience :: Science/Research',
'License :: OSI Approved :: BSD License',
'Topic :: Scientific/Engineering']
)
vendor/oosuite/OpenOpt/openopt.egg-info/dependency_links.txt 000644 000766 000024 00000000001 12100432771 025430 0 ustar 00rfadams staff 000000 000000
vendor/oosuite/OpenOpt/openopt.egg-info/not-zip-safe 000644 000766 000024 00000000001 12100432771 023610 0 ustar 00rfadams staff 000000 000000
vendor/oosuite/OpenOpt/openopt.egg-info/PKG-INFO 000644 000766 000024 00000001122 12100432771 022453 0 ustar 00rfadams staff 000000 000000 Metadata-Version: 1.0
Name: openopt
Version: 0.43
Summary: A python module for numerical optimization
Home-page: ('http://openopt.org',)
Author: ('Dmitrey Kroshko',)
Author-email: ('dmitrey-at-openopt-dot-org',)
License: new BSD
Download-URL: http://openopt.org/images/3/33/OpenOpt.zip
Description:
Platform: UNKNOWN
Classifier: Development Status :: 4 - Beta
Classifier: Environment :: Console
Classifier: Intended Audience :: Developers
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: BSD License
Classifier: Topic :: Scientific/Engineering
vendor/oosuite/OpenOpt/openopt.egg-info/requires.txt 000644 000766 000024 00000000022 12100432771 023754 0 ustar 00rfadams staff 000000 000000 numpy
setproctitle vendor/oosuite/OpenOpt/openopt.egg-info/SOURCES.txt 000644 000766 000024 00000016022 12100432771 023247 0 ustar 00rfadams staff 000000 000000 README.txt
setup.cfg
setup.py
openopt/__init__.py
openopt/oo.py
openopt/ooVersionNumber.py
openopt.egg-info/PKG-INFO
openopt.egg-info/SOURCES.txt
openopt.egg-info/dependency_links.txt
openopt.egg-info/not-zip-safe
openopt.egg-info/requires.txt
openopt.egg-info/top_level.txt
openopt/examples/GUI_1.py
openopt/examples/__init__.py
openopt/examples/dfp_1.py
openopt/examples/dfp_2.py
openopt/examples/eig_1.py
openopt/examples/eig_2.py
openopt/examples/glp_1.py
openopt/examples/glp_2.py
openopt/examples/glp_3.py
openopt/examples/glp_Ab_c.py
openopt/examples/lcp_1.py
openopt/examples/llavp_1.py
openopt/examples/llsp_1.py
openopt/examples/llsp_2.py
openopt/examples/lp_1.py
openopt/examples/lunp_1.py
openopt/examples/mcp_1.py
openopt/examples/milp_1.py
openopt/examples/minlp_1.py
openopt/examples/miqcqp_1.py
openopt/examples/miqp_1.py
openopt/examples/mmp_1.py
openopt/examples/mmp_2.py
openopt/examples/nllsp_1.py
openopt/examples/nlp_1.py
openopt/examples/nlp_11.py
openopt/examples/nlp_2.py
openopt/examples/nlp_3.py
openopt/examples/nlp_4.py
openopt/examples/nlp_5.py
openopt/examples/nlp_ALGENCAN.py
openopt/examples/nlp_bench_1.py
openopt/examples/nlp_bench_2.py
openopt/examples/nlp_d2f.py
openopt/examples/nlsp_1.py
openopt/examples/nlsp_constrained.py
openopt/examples/nsp_1.py
openopt/examples/nssolveVSfsolve_1.py
openopt/examples/qcqp_1.py
openopt/examples/qp_1.py
openopt/examples/sdp_1.py
openopt/examples/sle_1.py
openopt/examples/snle_1.py
openopt/examples/socp_1.py
openopt/examples/stab_1.py
openopt/examples/tsp_1.py
openopt/examples/tsp_2.py
openopt/examples/tsp_3.py
openopt/examples/tsp_4.py
openopt/kernel/DFP.py
openopt/kernel/EIG.py
openopt/kernel/GLP.py
openopt/kernel/GUI.py
openopt/kernel/IP.py
openopt/kernel/LCP.py
openopt/kernel/LLAVP.py
openopt/kernel/LLSP.py
openopt/kernel/LP.py
openopt/kernel/LUNP.py
openopt/kernel/MCP.py
openopt/kernel/MILP.py
openopt/kernel/MINLP.py
openopt/kernel/MMP.py
openopt/kernel/MOP.py
openopt/kernel/NLLSP.py
openopt/kernel/NLP.py
openopt/kernel/NLSP.py
openopt/kernel/NSP.py
openopt/kernel/ODE.py
openopt/kernel/Point.py
openopt/kernel/QP.py
openopt/kernel/SDP.py
openopt/kernel/SLE.py
openopt/kernel/SOCP.py
openopt/kernel/STAB.py
openopt/kernel/TSP.py
openopt/kernel/__init__.py
openopt/kernel/baseProblem.py
openopt/kernel/baseSolver.py
openopt/kernel/check.py
openopt/kernel/fdmisc.py
openopt/kernel/iterPrint.py
openopt/kernel/mfa.py
openopt/kernel/nonLinFuncs.py
openopt/kernel/nonOptMisc.py
openopt/kernel/ooIter.py
openopt/kernel/ooMisc.py
openopt/kernel/ooRun.py
openopt/kernel/oographics.py
openopt/kernel/oologfcn.py
openopt/kernel/pointProjection.py
openopt/kernel/residuals.py
openopt/kernel/runProbSolver.py
openopt/kernel/setDefaultIterFuncs.py
openopt/solvers/__init__.py
openopt/solvers/BrasilOpt/__init__.py
openopt/solvers/BrasilOpt/algencan_oo.py
openopt/solvers/CVXOPT/CVXOPT_LLSP_Solver.py
openopt/solvers/CVXOPT/CVXOPT_LP_Solver.py
openopt/solvers/CVXOPT/CVXOPT_QP_Solver.py
openopt/solvers/CVXOPT/CVXOPT_SDP_Solver.py
openopt/solvers/CVXOPT/CVXOPT_SOCP_Solver.py
openopt/solvers/CVXOPT/__init__.py
openopt/solvers/CVXOPT/cvxopt_lp_oo.py
openopt/solvers/CVXOPT/cvxopt_misc.py
openopt/solvers/CVXOPT/cvxopt_qp_oo.py
openopt/solvers/CVXOPT/cvxopt_sdp_oo.py
openopt/solvers/CVXOPT/cvxopt_socp_oo.py
openopt/solvers/CVXOPT/dsdp_oo.py
openopt/solvers/CVXOPT/glpk_oo.py
openopt/solvers/CoinOr/__init__.py
openopt/solvers/CoinOr/cplex_oo.py
openopt/solvers/CoinOr/ipopt_oo.py
openopt/solvers/HongKongOpt/LCPSolve.py
openopt/solvers/HongKongOpt/QPSolve.py
openopt/solvers/HongKongOpt/QPSparse.py
openopt/solvers/HongKongOpt/__init__.py
openopt/solvers/HongKongOpt/lcpsolve_oo.py
openopt/solvers/HongKongOpt/qlcp.py
openopt/solvers/HongKongOpt/qlcp_oo.py
openopt/solvers/HongKongOpt/sqlcp.py
openopt/solvers/HongKongOpt/sqlcp_oo.py
openopt/solvers/Standalone/__init__.py
openopt/solvers/Standalone/arpack_oo.py
openopt/solvers/Standalone/asa_oo.py
openopt/solvers/Standalone/branb_oo.py
openopt/solvers/Standalone/bvls_oo.py
openopt/solvers/Standalone/defaultSLEsolver_oo.py
openopt/solvers/Standalone/galileo_oo.py
openopt/solvers/Standalone/lsqr.py
openopt/solvers/Standalone/lsqr_oo.py
openopt/solvers/Standalone/pclp_oo.py
openopt/solvers/Standalone/pointProjection.py
openopt/solvers/Standalone/pswarm_oo.py
openopt/solvers/Standalone/pymls_oo.py
openopt/solvers/Standalone/sa.py
openopt/solvers/Standalone/sa_oo.py
openopt/solvers/Standalone/tsp.py
openopt/solvers/UkrOpt/Dilation.py
openopt/solvers/UkrOpt/Dilation2.py
openopt/solvers/UkrOpt/Dilation3.py
openopt/solvers/UkrOpt/PolytopProjection.py
openopt/solvers/UkrOpt/ShorEllipsoid_oo.py
openopt/solvers/UkrOpt/UkrOptMisc.py
openopt/solvers/UkrOpt/__init__.py
openopt/solvers/UkrOpt/amsg2p.py
openopt/solvers/UkrOpt/amsg2p_oo.py
openopt/solvers/UkrOpt/de_oo.py
openopt/solvers/UkrOpt/goldenSection_oo.py
openopt/solvers/UkrOpt/gsubg_oo.py
openopt/solvers/UkrOpt/ii_engine.py
openopt/solvers/UkrOpt/interalgCons.py
openopt/solvers/UkrOpt/interalgLLR.py
openopt/solvers/UkrOpt/interalgMOP.py
openopt/solvers/UkrOpt/interalgMisc.py
openopt/solvers/UkrOpt/interalgODE.py
openopt/solvers/UkrOpt/interalgT.py
openopt/solvers/UkrOpt/interalg_oo.py
openopt/solvers/UkrOpt/lincher_oo.py
openopt/solvers/UkrOpt/nsmm_oo.py
openopt/solvers/UkrOpt/nssolve_oo.py
openopt/solvers/UkrOpt/ralg_oo.py
openopt/solvers/lp_solve/__init__.py
openopt/solvers/lp_solve/lpSolve_oo.py
openopt/solvers/nlopt/NLOPT_AUX.py
openopt/solvers/nlopt/NLOPT_BASE.py
openopt/solvers/nlopt/__init__.py
openopt/solvers/nlopt/auglag_oo.py
openopt/solvers/nlopt/bobyqa_oo.py
openopt/solvers/nlopt/isres_oo.py
openopt/solvers/nlopt/mlsl_oo.py
openopt/solvers/nlopt/mma_oo.py
openopt/solvers/nlopt/ptn_oo.py
openopt/solvers/nlopt/sbplx_oo.py
openopt/solvers/nlopt/slmvm1_oo.py
openopt/solvers/nlopt/slmvm2_oo.py
openopt/solvers/nlopt/stogo_oo.py
openopt/solvers/scipy_optim/__init__.py
openopt/solvers/scipy_optim/lapack_dgelss_oo.py
openopt/solvers/scipy_optim/lapack_sgelss_oo.py
openopt/solvers/scipy_optim/lsmr_oo.py
openopt/solvers/scipy_optim/numpy_eig_oo.py
openopt/solvers/scipy_optim/scipy_anderson2_oo.py
openopt/solvers/scipy_optim/scipy_anderson_oo.py
openopt/solvers/scipy_optim/scipy_bfgs_oo.py
openopt/solvers/scipy_optim/scipy_broyden1_oo.py
openopt/solvers/scipy_optim/scipy_broyden2_oo.py
openopt/solvers/scipy_optim/scipy_broyden3_oo.py
openopt/solvers/scipy_optim/scipy_broyden_generalized_oo.py
openopt/solvers/scipy_optim/scipy_cg_oo.py
openopt/solvers/scipy_optim/scipy_cobyla_oo.py
openopt/solvers/scipy_optim/scipy_fmin_oo.py
openopt/solvers/scipy_optim/scipy_fminbound_oo.py
openopt/solvers/scipy_optim/scipy_fsolve_oo.py
openopt/solvers/scipy_optim/scipy_lbfgsb_oo.py
openopt/solvers/scipy_optim/scipy_leastsq_oo.py
openopt/solvers/scipy_optim/scipy_ncg_oo.py
openopt/solvers/scipy_optim/scipy_powell_oo.py
openopt/solvers/scipy_optim/scipy_slsqp_oo.py
openopt/solvers/scipy_optim/scipy_tnc_oo.py
openopt/tests/__init__.py
openopt/tests/blockMisc.py
openopt/tests/chain.py
openopt/tests/nlp1.py
openopt/tests/nlp2.py
openopt/tests/nlpLC.py
openopt/tests/nlsp1.py
openopt/tests/nsp1.py
openopt/tests/nsp2.py
openopt/tests/oosolver.py
openopt/tests/rastrigin.py vendor/oosuite/OpenOpt/openopt.egg-info/top_level.txt 000644 000766 000024 00000000010 12100432771 024103 0 ustar 00rfadams staff 000000 000000 openopt
vendor/oosuite/OpenOpt/openopt/.___init__.py 000755 000766 000024 00000000324 12100431107 022210 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/openopt/__init__.py 000755 000766 000024 00000002137 12100431107 021777 0 ustar 00rfadams staff 000000 000000 #! /usr/bin/env python
#from .ooVersionNumber import __version__
import os, sys
curr_dir = ''.join([elem + os.sep for elem in __file__.split(os.sep)[:-1]])
sys.path += [curr_dir, curr_dir + 'kernel']
from ooVersionNumber import __version__
from oo import *
#from kernel.GUI import manage
#from kernel.oologfcn import OpenOptException
#from kernel.nonOptMisc import oosolver
from GUI import manage
from oologfcn import OpenOptException
from nonOptMisc import oosolver
from mfa import MFA
isE = False
try:
import enthought
isE = True
except ImportError:
pass
try:
import envisage
import mayavi
isE = True
except ImportError:
pass
try:
import xy
isE = False
except ImportError:
pass
if isE:
s = """
Seems like you are using OpenOpt from
commercial Enthought Python Distribution;
consider using free GPL-licensed alternatives
PythonXY (http://www.pythonxy.com) or
Sage (http://sagemath.org) instead.
"""
print(s)
#__all__ = filter(lambda s:not s.startswith('_'),dir())
#from numpy.testing import NumpyTest
#test = NumpyTest().test
vendor/oosuite/OpenOpt/openopt/._doc 000755 000766 000024 00000000324 12100431107 020567 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/openopt/doc/ 000755 000766 000024 00000000000 12100431107 020425 5 ustar 00rfadams staff 000000 000000 vendor/oosuite/OpenOpt/openopt/._examples 000755 000766 000024 00000000324 12100431107 021640 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/openopt/examples/ 000755 000766 000024 00000000000 12100431107 021476 5 ustar 00rfadams staff 000000 000000 vendor/oosuite/OpenOpt/openopt/._kernel 000755 000766 000024 00000000324 12100431107 021302 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/openopt/kernel/ 000755 000766 000024 00000000000 12100431107 021140 5 ustar 00rfadams staff 000000 000000 vendor/oosuite/OpenOpt/openopt/._oo.py 000755 000766 000024 00000000324 12100431107 021066 0 ustar 00rfadams staff 000000 000000 Mac OS X 2 ATTR < < com.apple.quarantine q/0002;5102310c;Safari;E78377CF-5ECF-413D-BFCA-810632705AF0 vendor/oosuite/OpenOpt/openopt/oo.py 000755 000766 000024 00000070252 12100431107 020660 0 ustar 00rfadams staff 000000 000000 import os,sys
sys.path.append(os.getcwd()+os.sep+'kernel')
from LP import LP as CLP
from LCP import LCP as CLCP
from EIG import EIG as CEIG
from SDP import SDP as CSDP
from QP import QP as CQP
from MILP import MILP as CMILP
from STAB import STAB as CSTAB
from MCP import MCP as CMCP
from TSP import TSP as CTSP
from NSP import NSP as CNSP
from NLP import NLP as CNLP
from MOP import MOP as CMOP
from MINLP import MINLP as CMINLP
from NLSP import NLSP as CNLSP
from NLLSP import NLLSP as CNLLSP
from GLP import GLP as CGLP
from SLE import SLE as CSLE
from LLSP import LLSP as CLLSP
from MMP import MMP as CMMP
from LLAVP import LLAVP as CLLAVP
from LUNP import LUNP as CLUNP
from SOCP import SOCP as CSOCP
from DFP import DFP as CDFP
from IP import IP as CIP
from ODE import ODE as CODE
def MILP(*args, **kwargs):
"""
MILP: constructor for Mixed Integer Linear Problem assignment
f' x -> min
subjected to
lb <= x <= ub
A x <= b
Aeq x = beq
for all i from intVars: i-th coordinate of x is required to be integer
for all j from binVars: j-th coordinate of x is required to be from {0, 1}
Examples of valid calls:
p = MILP(f, )
p = MILP(f=objFunVector, )
p = MILP(f, A=A, intVars = myIntVars, Aeq=Aeq, b=b, beq=beq, lb=lb, ub=ub, binVars = binVars)
See also: /examples/milp_*.py
:Parameters:
- intVars : Python list of those coordinates that are required to be integers.
- binVars : Python list of those coordinates that are required to be binary.
all other input parameters are same to LP class constructor ones
:Returns:
OpenOpt MILP class instance
Notes
-----
Solving of MILPs is performed via
r = p.solve(string_name_of_solver)
or p.maximize, p.minimize
r.xf - desired solution (NaNs if a problem occured)
r.ff - objFun value () (NaN if a problem occured)
(see also other r fields)
Solvers available for now:
lpSolve (LGPL) - requires lpsolve + Python bindings installations (all mentioned is available in http://sourceforge.net/projects/lpsolve)
glpk (GPL 2) - requires glpk + CVXOPT v >= 1.0 installations (read OO MILP webpage for more details)
"""
return CMILP(*args, **kwargs)
def LP(*args, **kwargs):
"""
LP: constructor for Linear Problem assignment
f' x -> min
subjected to
lb <= x <= ub
A x <= b
Aeq x = beq
valid calls are:
p = LP(f, )
p = LP(f=objFunVector, )
p = LP(f, A=A, Aeq=Aeq, Awhole=Awhole, b=b, beq=beq, bwhole=bwhole, dwhole=dwhole, lb=lb, ub=ub)
See also: /examples/lp_*.py
:Parameters:
f: vector of length n
A: size m1 x n matrix, subjected to A * x <= b
Aeq: size m2 x n matrix, subjected to Aeq * x = beq
b, beq: corresponding vectors of lengthes m1, m2
lb, ub: vectors of length n, some coords may be +/- inf
:Returns:
OpenOpt LP class instance
Notes
-----
Solving of LPs is performed via
r = p.solve(string_name_of_solver)
or p.maximize, p.minimize
r.xf - desired solution (NaNs if a problem occured)
r.ff - objFun value () (NaN if a problem occured)
(see also other r fields)
Solvers available for now:
pclp (BSD) - premature but pure Python implementation with permissive license
lpSolve (LGPL) - requires lpsolve + Python bindings installations (all mentioned is available in http://sourceforge.net/projects/lpsolve)
cvxopt_lp (GPL) - requires CVXOPT (http://abel.ee.ucla.edu/cvxopt)
glpk(GPL2) - requires CVXOPT(http://abel.ee.ucla.edu/cvxopt) & glpk (www.gnu.org/software/glpk)
converter to NLP. Example: r = p.solve('nlp:ipopt')
"""
return CLP(*args, **kwargs)
def STAB(*args, **kwargs):
"""
STAB: constructor for graph stability number problem assignment
valid calls are:
p = STAB(graph, )
p = STAB(graph = graph, )
See also: /examples/stab_*.py
:Parameters:
graph: networkx graph instance
:Returns:
OpenOpt STAB class instance
Notes
-----
Solving of OpenOpt STAB problems is performed via
r = p.solve(string_name_of_solver)
or r = p.manage(string_name_of_solver) (to enable basic GUI)
r.solution - desired solution (python list of nodes)
r.ff - objFun value (NaN if a problem occured)
(see also other r fields)
Solvers available for now:
* interalg (BSD license)
* MILP solvers at http://openopt.org/MILP
"""
return CSTAB(*args, **kwargs)
def MCP(*args, **kwargs):
"""
MCP: constructor for maximum clique problem assignment
valid calls are:
p = MCP(graph, )
p = MCP(graph = graph, )
See also: /examples/mcp_*.py
:Parameters:
graph: networkx graph instance
:Returns:
OpenOpt MCP class instance
Notes
-----
Solving of OpenOpt MCP problems is performed via
r = p.solve(string_name_of_solver)
or r = p.manage(string_name_of_solver) (to enable basic GUI)
r.solution - desired solution (python list of nodes)
r.ff - objFun value (NaN if a problem occured)
(see also other r fields)
Solvers available for now:
* interalg (BSD license)
* MILP solvers at http://openopt.org/MILP
"""
return CMCP(*args, **kwargs)
def TSP(*args, **kwargs):
"""
TSP: constructor for traveling salesman problem assignment
valid calls are:
p = TSP(graph, )
p = TSP(graph = graph, )
See also: /examples/tsp_*.py
:Parameters:
graph: networkx graph instance
:Returns:
OpenOpt TSP class instance
Notes
-----
Solving of OpenOpt TSP problems is performed via
r = p.solve(string_name_of_solver)
or r = p.manage(string_name_of_solver) (to enable basic GUI)
r.solution - desired solution (python list of edges)
r.ff - objFun value (NaN if a problem occured)
(see also other r fields)
Solvers available for now:
* interalg (BSD license)
* MILP solvers at http://openopt.org/MILP
"""
return CTSP(*args, **kwargs)
def LCP(*args, **kwargs):
"""
LCP: constructor for Linear Complementarity Problem assignment
find w, z: w = Mz + q
valid calls are:
p = LCP(M, q, )
p = LCP(M=M, q=q, )
See also: /examples/lcp_*.py
:Parameters:
M: numpy array of size n x n
q: vector of length n
:Returns:
OpenOpt LCP class instance
Notes
-----
Solving of LCPs is performed via
r = p.solve(string_name_of_solver)
r.xf - desired solution (1st n/2 coords are w, other n/2 coords are z)
r.ff - objFun value (max residual of Mz+q-w) (NaN if a problem occured)
(see also other r fields)
Solvers available for now:
lcpsolve (BSD)
"""
return CLCP(*args, **kwargs)
def EIG(*args, **kwargs):
"""
EIG: constructor for Eigenvalues Problem assignment
to solve standard eigenvalue problem:
find eigenvalues and eigenvectors of square matrix A:
A x = lambda x
or general eigenvalue problem:
A x = lambda M x
valid calls are:
p = EIG(M, q, )
p = EIG(M=M, q=q, )
See also: /examples/eig_*.py
:Parameters:
A, (optional) M: numpy array or scipy sparse matrix of size n x n
:Returns:
OpenOpt EIG class instance
Notes
-----
Solving of EIGs is performed via
r = p.solve(string_name_of_solver)
see http://openopt.org/EIG for more info
Solvers available for now:
arpack (license: BSD)
"""
return CEIG(*args, **kwargs)
def SDP(*args, **kwargs):
"""
SDP: constructor for SemiDefinite Problem assignment
f' x -> min
subjected to
lb <= x <= ub
A x <= b
Aeq x = beq
For all i = 0, ..., I: Sum [j = 0, ..., n-1] {S_i_j x_j} <= d_i (matrix componentwise inequality),
d_i are square matrices
S_i_j are square positive semidefinite matrices of size same to d_i
valid calls are:
p = SDP(f, )
p = SDP(f=objFunVector, )
p = SDP(f, S=S, d=d, A=A, Aeq=Aeq, Awhole=Awhole, b=b, beq=beq, bwhole=bwhole, dwhole=dwhole, lb=lb, ub=ub)
See also: /examples/sdp_*.py
:Parameters:
f: vector of length n
S: Python dict of square matrices S[0, 0], S[0,1], ..., S[I,J]
S[i, j] are real symmetric positive-definite matrices
d: Python dict of square matrices d[0], ..., d[I]
A: size m1 x n matrix, subjected to A * x <= b
Aeq: size m2 x n matrix, subjected to Aeq * x = beq
b, beq: corresponding vectors of lengthes m1, m2
lb, ub: vectors of length n, some coords may be +/- inf
:Returns:
OpenOpt SDP class instance
Notes
-----
Solving of SDPs is performed via
r = p.solve(string_name_of_solver)
r.xf - desired solution (NaNs if a problem occured)
r.ff - objFun value () (NaN if a problem occured)
(see also other r fields)
Solvers available for now:
cvxopt_sdp (LGPL) - requires CVXOPT (http://abel.ee.ucla.edu/cvxopt)
dsdp (GPL) - requires CVXOPT + DSDP installation, can't handle linear equality constraints Aeq x = beq
"""
return CSDP(*args, **kwargs)
def SOCP(*args, **kwargs):
"""
SOCP: constructor for Second-Order Cone Problem assignment
f' x -> min
subjected to
lb <= x <= ub
Aeq x = beq
For all i = 0, ..., I: ||C_i x + d_i|| <= q_i x + s_i
valid calls are:
p = SDP(f, )
p = SDP(f=objFunVector, )
p = SDP(f, S=S, d=d, A=A, Aeq=Aeq, Awhole=Awhole, b=b, beq=beq, bwhole=bwhole, dwhole=dwhole, lb=lb, ub=ub)
See also: /examples/sdp_*.py
:Parameters:
f: vector of length n
Aeq: size M x n matrix, subjected to Aeq * x = beq
beq: corresponding vector of length M
C: Python list of matrices C_i of shape (m_i, n)
d: Python list of vectors of length m_i
q: Python list of vectors of length n
s: Python list of numbers, len(s) = n
lb, ub: vectors of length n, some coords may be +/- inf
:Returns:
OpenOpt SDP class instance
Notes
-----
Solving of SOCPs is performed via
r = p.solve(string_name_of_solver)
r.xf - desired solution (NaNs if a problem occured)
r.ff - objFun value () (NaN if a problem occured)
(see also other r fields)
Solvers available for now:
cvxopt_socp (LGPL) - requires CVXOPT (http://abel.ee.ucla.edu/cvxopt)
"""
return CSOCP(*args, **kwargs)
def QP(*args, **kwargs):
"""
QP: constructor for Quadratic Problem assignment
1/2 x' H x + f' x -> min
subjected to
A x <= b
Aeq x = beq
lb <= x <= ub
Examples of valid calls:
p = QP(H, f, )
p = QP(numpy.ones((3,3)), f=numpy.array([1,2,4]), )
p = QP(f=range(8)+15, H = numpy.diag(numpy.ones(8)), )
p = QP(H, f, A=A, Aeq=Aeq, b=b, beq=beq, lb=lb, ub=ub, )
See also: /examples/qp_*.py
INPUT:
H: size n x n matrix, symmetric, positive-definite
f: vector of length n
lb, ub: vectors of length n, some coords may be +/- inf
A: size m1 x n matrix, subjected to A * x <= b
Aeq: size m2 x n matrix, subjected to Aeq * x = beq
b, beq: vectors of lengths m1, m2
Alternatively to A/Aeq you can use Awhole matrix as it's described in LP documentation (or both A, Aeq, Awhole)
OUTPUT: OpenOpt QP class instance
Solving of QPs is performed via
r = p.solve(string_name_of_solver)
r.xf - desired solution (NaNs if a problem occured)
r.ff - objFun value (NaN if a problem occured)
(see also other r fields)
Solvers available for now:
cvxopt_qp (GPL) - requires CVXOPT (http://abel.ee.ucla.edu/cvxopt)
converter to NLP. Example: r = p.solve('nlp:ipopt')
"""
return CQP(*args, **kwargs)
def NLP(*args, **kwargs):
"""
NLP: constructor for general Non-Linear Problem assignment
f(x) -> min (or -> max)
subjected to
c(x) <= 0
h(x) = 0
A x <= b
Aeq x = beq
lb <= x <= ub
Examples of valid usage:
p = NLP(f, x0, )
p = NLP(f=objFun, x0 = myX0, )
p = NLP(f, x0, A=A, df = objFunGradient, Aeq=Aeq, b=b, beq=beq, lb=lb, ub=ub)
See also: /examples/nlp_*.py
INPUTS:
f: objFun
x0: start point, vector of length n
Optional:
name: problem name (string), is used in text & graphics output
df: user-supplied gradient of objective function
c, h - functions defining nonlinear equality/inequality constraints
dc, dh - functions defining 1st derivatives of non-linear constraints
A: size m1 x n matrix, subjected to A * x <= b
Aeq: size m2 x n matrix, subjected to Aeq * x = beq
b, beq: corresponding vectors of lengthes m1, m2
lb, ub: vectors of length n subjected to lb <= x <= ub constraints, may include +/- inf values
iprint = {10}: print text output every iteration
goal = {'minimum'} | 'min' | 'maximum' | 'max' - minimize or maximize objective function
diffInt = {1e-7} : finite-difference gradient aproximation step, scalar or vector of length nVars
scale = {None} : scale factor, see /examples/badlyScaled.py for more details
stencil = {1}|2|3: finite-differences derivatives approximation stencil,
used by most of solvers (except scipy_cobyla) when no user-supplied for
objfun / nonline constraints derivatives are provided
1: (f(x+dx)-f(x))/dx (faster but less precize)
2: (f(x+dx)-f(x-dx))/(2*dx) (slower but more exact)
3: (-f(x+2*dx)+8*f(x+dx)-8*f(x-dx)+f(x-2*dx))/(12*dx) (even more slower, but even more exact)
check.df, check.dc, check.dh: if set to True, OpenOpt will check user-supplied gradients.
args (or args.f, args.c, args.h) - additional arguments to objFunc and non-linear constraints,
see /examples/userArgs.py for more details.
contol: max allowed residual in optim point
(for any constraint from problem constraints:
constraint(x_optim) < contol is required from solver)
stop criteria:
maxIter {400}
maxFunEvals {1e5}
maxCPUTime {inf}
maxTime {inf}
maxLineSearch {500}
fEnough {-inf for min problems, +inf for max problems}:
stop if objFunc vulue better than fEnough and all constraints less than contol
ftol {1e-6}: used in stop criterium || f[iter_k] - f[iter_k+1] || < ftol
xtol {1e-6}: used in stop criterium || x[iter_k] - x[iter_k+1] || < xtol
gtol {1e-6}: used in stop criteria || gradient(x[iter_k]) || < gtol
callback - user-defined callback function(s), see /examples/userCallback.py
Notes:
1) for more safety default values checking/reassigning (via print p.maxIter / prob.maxIter = 400) is recommended
(they may change in future OpenOpt versions and/or not updated in time in the documentation)
2) some solvers may ignore some of the stop criteria above and/or use their own ones
3) for NSP constructor ftol, xtol, gtol defaults may have other values
graphic options:
plot = {False} | True : plot figure (now implemented for UC problems only), requires matplotlib installed
color = {'blue'} | black | ... (any valid matplotlib color)
specifier = {'-'} | '--' | ':' | '-.' - plot specifier
show = {True} | False : call pylab.show() after solver finish or not
xlim {(nan, nan)}, ylim {(nan, nan)} - initial estimation for graphical output borders
(you can use for example p.xlim = (nan, 10) or p.ylim = [-8, 15] or p.xlim=[inf, 15], only real finite values will be taken into account)
for constrained problems ylim affects only 1st subplot
p.graphics.xlabel or p.xlabel = {'time'} | 'cputime' | 'iter' # desired graphic output units in x-axe, case-unsensetive
Note: some Python IDEs have problems with matplotlib!
Also, after assignment NLP instance you may modify prob fields inplace:
p.maxIter = 1000
p.df = lambda x: cos(x)
OUTPUT: OpenOpt NLP class instance
Solving of NLPs is performed via
r = p.solve(string_name_of_solver)
or p.maximize, p.minimize
r.xf - desired solution (NaNs if a problem occured)
r.ff - objFun value (NaN if a problem occured)
(see also other fields, such as CPUTimeElapsed, TimeElapsed, isFeasible, iter etc, via dir(r))
Solvers available for now:
single-variable:
goldenSection, scipy_fminbound (latter is not recommended)
(both these solvers require finite lb-ub and ignore user-supplied gradient)
unconstrained:
scipy_bfgs, scipy_cg, scipy_ncg,
(these ones cannot handle user-provided gradient) scipy_powell and scipy_fmin
amsg2p - requires knowing fOpt (optimal value)
box-bounded:
scipy_lbfgsb, scipy_tnc - require scipy installed
bobyqa - doesn't use derivatives; requires http://openopt.org/nlopt installed
ptn, slmvm1, slmvm2 - require http://openopt.org/nlopt installed
all constraints:
ralg
ipopt (requires ipopt + pyipopt installed)
scipy_slsqp
scipy_cobyla (this one cannot handle user-supplied gradients)
lincher (requires CVXOPT QP solver),
gsubg - for large-scaled problems
algencan (ver. 2.0.3 or more recent, very powerful constrained solver, GPL,
requires ALGENCAN + Python interface installed,
see http://www.ime.usp.br/~egbirgin/tango/)
mma and auglag - require http://openopt.org/nlopt installed
"""
return CNLP(*args, **kwargs)
def MINLP(*args, **kwargs):
"""
MINLP: constructor for general Mixed-Integer Non-Linear Problem assignment
parameters and usage: same to NLP, + parameters
discreteVars: dictionary numberOfCoord <-> list (or tuple) of allowed values, eg
p.discreteVars = {0: [1, 2.5], 15: (3.1, 4), 150: [4,5, 6]}
discrtol (default 1e-5) - tolerance required for discrete constraints
available solvers:
branb (branch-and-bound) - translation of fminconset routine, requires non-default string parameter nlpSolver
"""
return CMINLP(*args, **kwargs)
def NSP(*args, **kwargs):
"""
Non-Smooth Problem constructor
Same usage as NLP (see help(NLP) and /examples/nsp_*.py), but default values of contol, xtol, ftol, diffInt may differ
Also, default finite-differences derivatives approximation stencil is 3 instead of 1 for NLP
Solvers available for now:
ralg - all constraints, medium-scaled (nVars = 1...1000), can handle user-provided gradient/subgradient
amsg2p - requires knowing fOpt (optimal value), medium-scaled (nVars = 1...1000), can handle user-provided gradient/subgradient
gsubg - for large-scaled problems
scipy_fmin - a Nelder-Mead simplex algorithm implementation, cannot handle constraints and derivatives
sbplx - A variant of Nelder-Mead algorithm; requires http://openopt.org/nlopt installed
ShorEllipsoid (unconstrained for now) - small-scale, nVars=1...10, requires r0: ||x0-x*||<=r0
"""
return CNSP(*args, **kwargs)
def NLSP(*args, **kwargs):
"""
Solving systems of n non-linear equations with n variables
Parameters and usage: same as NLP
(see help(NLP) and /examples/nlsp_*.py)
Solvers available for now:
scipy_fsolve (can handle df);
converter to NLP. Example: r = p.solve('nlp:ipopt');
nssolve (primarily for non-smooth and noisy funcs; can handle all types of constraints and 1st derivatives df,dc,dh; splitting equations to Python list or tuple is recommended to speedup calculations)
(these ones below are very unstable and can't use user-supplied gradient - at least, for scipy 0.6.0)
scipy_anderson
scipy_anderson2
scipy_broyden1
scipy_broyden2
scipy_broyden3
scipy_broyden_generalized
"""
r = CNLSP(*args, **kwargs)
r.pWarn('''
OpenOpt NLSP class had been renamed to SNLE
(system of nonlinear equations), use "SNLE" instead of "NLSP"
''')
return r
def SNLE(*args, **kwargs):
"""
Solving systems of n non-linear equations with n variables
Parameters and usage: same as NLP
(see help(NLP) and /examples/nlsp_*.py)
Solvers available for now:
scipy_fsolve (can handle df);
converter to NLP. Example: r = p.solve('nlp:ipopt');
nssolve (primarily for non-smooth and noisy funcs; can handle all types of constraints and 1st derivatives df,dc,dh; splitting equations to Python list or tuple is recommended to speedup calculations)
(these ones below are very unstable and can't use user-supplied gradient - at least, for scipy 0.6.0)
scipy_anderson
scipy_anderson2
scipy_broyden1
scipy_broyden2
scipy_broyden3
scipy_broyden_generalized
"""
return CNLSP(*args, **kwargs)
def NLLSP(*args, **kwargs):
"""
Given set of non-linear equations
f1(x)=0, f2(x)=0, ... fm(x)=0
search for x: f1(x, )^2 + ,,, + fm(x, )^2 -> min
Parameters and usage: same as NLP
(see help(openopt.NLP) and /examples/nllsp_*.py)
Solvers available for now:
scipy_leastsq (requires scipy installed)
converter to NLP. Example: r = p.solve('nlp:ralg')
"""
return CNLLSP(*args, **kwargs)
def MOP(*args, **kwargs):
'''
Multiobjective optimization
Search for weak or strong Pareto front
Solvers available for now:
interalg (http://openopt.org/interalg)
'''
return CMOP(*args, **kwargs)
def IP(*args, **kwargs):
"""
Integrate a function f: R^n -> R over a given domain lb_i <= x_i <= ub_i
"""
return CIP(*args, **kwargs)
def ODE(*args, **kwargs):
"""
Solve ODE dy/dt = f(y,t), y(0) = y0
"""
return CODE(*args, **kwargs)
def SLE(*args, **kwargs):
"""
SLE: constructor for system of linear equations C*x = d assignment
Examples of valid usage:
p = SLE(C, d,